NAb
NAE
DM
CARNEQIE:
PITTSBURGH, PA 15213-3690
Computational Modeling
and Mathematics
Applied to
the Physical Sciences
Committee on the Applications of Mathematics
Office of Mathematical Sciences
Commission on Physical Sciences, Mathematics, and Resources
National Research Council
NATIONAL ACADEMY PRESS
Washington, D.C. 1984
NOTICE: The project that is the subject of this report was approved by
the Governing Board of the National Research Council, whose members
are drawn from the councils of the National Academy of Sciences, the
National Academy of Engineering, and the Institute of Medicine. The
members of the committee responsible for the report were chosen for their
special competences and with regard for appropriate balance.
This report has been reviewed by a group other than the authors ac-
cording to procedures approved by a Report Review Committee consisting
of members of the National Academy of Sciences, the National Academy
of Engineering, and the Institute of Medicine.
The National Research Council was established by the National
Academy of Sciences in 1916 to associate the broad community of science
and technology with the Academy's purposes of furthering knowledge and
of advising the federal government. The Council operates in accordance
with general policies determined by the Academy under the authority of
its congressional charter of 1863, which establishes the Academy as a
private, nonprofit, self-governing membership corporation. The Council
has become the principal operating agency of both the National Academy
of Sciences and the National Academy of Engineering in the conduct of
their sendees to the government, the public, and the scientific and en-
gineering communities. It is administered jointly by both Academies and
the Institute of Medicine. The National Academy of Engineering and the
Institute of Medicine were established in 1964 and 1970, respectively, un-
der the charter of the National Academy of Sciences.
This work was sponsored by the U.S. Department of Energy under
contract number DE FG1081ER12001.
Available from
Office of Mathematical Sciences
National Research Council
2101 Constitution Avenue, NW
Washington, DC 20418
COMMITTEE ON APPLICATIONS OF MATHEMATICS
WERNER C. RHEINBOLDT, University of Pittsburgh, Chairman
JAMES G. GLIMM, New York University
JAMES M. HYMAN, Los Alamos National Laboratory
ROBERT J. KEE, JR., Sandia National Laboratories
JAMES A. KRUMHANSL, Cornell University
JOHN E. OSBORN, University of Maryland
MARTIN H. SCHULTZ, Yale University
IVAR STAKGOLD, University of Delaware
JACOB K. GOLDHABER, Executive Secretary
The Committee on Applications of Mathematics is a standing committee of
the Commission on Physical Sciences, Mathematics, and Resources of the
National Research Council. It is charged with taking all actions deemed
appropriate to enhance the effectiveness of mathematics and mathematical
statistics in their applications.
m
i-n*if%f
COMMISSION ON PHYSICAL SCIENCES, MATHEMATICS, AND RESOURCI
HERBERT FRIEDMAN, National Research Council, Chairman
ELKAN R. BLOUT, Harvard Medical School
WILLIAM BROWDER, Princeton University
BERNARD F. BURKE, Massachusetts Institute of Technology
HERMAN CHERNOFF, Massachusetts Institute of Technology
MILDRED S. DRESSELHAUS, Massachusetts Institute of Technology
WALTER R. ECKELMANN, Sohio Petroleum Company
JOSEPH L. FISHER, Office of the Governor, Commonwealth of Virginia
JAMES C. FLETCHER, University of Pittsburgh
WILLIAM A. FOWLER, California Institute of Technology
GERHART FRIEDLANDER, Brookhaven National Laboratory
EDWARD A. FRIEMAN, Science Applications, Inc.
EDWARD D. GOLDBERG, Scripps Institution of Oceanography
CHARLES L. HOSLER, JR., Pennsylvania State University
KONRAD B. KRAUSKOPF, Stanford University
CHARLES J. MANKIN, Oklahoma Geological Survey
WALTER H. MUNK, University of California, San Diego
GEORGE E. PAKE, Xerox Research Center
ROBERT E. SIEVERS, University of Colorado
HOWARD E. SIMMONS, JR., E. I. du Pont de Nemours & Co., Inc.
JOHN D. SPENGLER, Harvard School of Public Health
HATTEN S. YODER, JR., Carnegie Institution of Washington
RAPHAEL G. KASPER, Executive Director
IV
CONTENTS
Overview 1
1. Introduction 3
2. Applications 8
2.1 Hydrodynamic Systems 8
2.2 Chemical Systems and Combustion 17
2.3 Plasma Physics 24
2.4 Particle Physics 26
2.5 Condensed-Matter Physics 27
2.6 Geophysical Applications 31
2.7 Meteorology 33
2.8 Astrophysical Applications 35
2.9 Structural Mechanics 37
2.10 Nondestructive Testing and Tomographic Reconstruction 38
2.11 Mathematical Models in the Biological Sciences 41
2.12 Electronic Components 45
3. Computational and Mathematical Difficulties 47
3.1 Degrees of Freedom 47
3.2 Different Time and Length Scales 49
3.3 Singularities in Coefficients, Data, or States 52
3.4 Boundary Conditions 52
3.5 Bifurcations and Chaos 54
3.6 El-Posed Inverse Problems 59
3.7 Effective Media 61
3.8 Validation, Error Assessments, and Sensitivity Analysis 63
4. Numerical Methods 67
4.1 Discretization Methods 67
4.2 Adaptive Grid Methods 72
4.3 Methods for Solving Discretized Differential Equations 75
4.4 Continuation and Homotopy Methods 78
4.5 Optimization Methods 80
4.6 Graph-Theoretical Methods 81
4.7 Splitting Methods and Defect Corrections 84
4.8 Monte Carlo Techniques 85
4.9 Problem-Dependent Methods 89
4.10 Computer Software and Hardware 93
Appendix A: Persons Contributing to this Report 96
OVERVIEW
The 1980s are a time of profound challenge to the technological strength
of the United States in the economic as well as the military spheres, and
our country's performance in research and development in its engineering
laboratories will be an important and perhaps determining aspect of our
success in meeting this challenge. Advanced engineering development is
now based mainly on scientific computation, which in turn relies on math-
ematical modeling and laboratory experiments. Together they represent
one of the areas in which the strength of nations is being tested today.
Mathematics is essential in the development of theoretical and com-
putational models for solving the highly complex problems of engineering
and basic science, which encompass a range of scientific difficulties. On one
side are questions of computer architecture and the science of algorithms.
On the other side is the modeling of chemical and physical processes
by means of mathematical equations. The issues are tied together by
mathematical theory, which seeks a full understanding of the nonlinear
phenomena contained in the equations and implements this understanding
through computational methods. This span of scientific activities forms
the subject that is known as applied mathematics.
The principal conclusion of this committee is that computational
modeling, which is a high-leverage element of our nation's scientific and
technological effort, requires increased emphasis and support. The conclu-
sion is documented by an examination of typical application areas, which
reveals the pervasive difficulties that accompany computation of realis-
tic problems and leads one to consider both what computers can do and
what they cannot currently do but might eventually be capable of doing.
As illustration, we examine several deep theoretical problems, including
turbulence and combustion. At the frontiers of attack on these prob-
lems we discover the limitations imposed by our current understanding of
model formulation and computational capability. We examine modeling
problems and algorithms adapted both to specific features of the desired
solution and to the computer architecture. We also examine computer
architecture and component design and manufacture as a mathematical
modeling problem.
The Committee recommends increased support for
1. Research in computational modeling and applied mathematics,
1
2. Computing facilities dedicated to this area,
3. Education and manpower development in computational and ap-
plied mathematics.
These recommended increases include financial support from govern-
ment and industry as well as institutional support from universities.
1. INTRODUCTION
The extensive use of computers in advanced development work began dur-
ing World War n, and today computing is a vital component of science,
engineering, and modern technology. Most advanced technological devel-
opment, from aircraft design to automobiles to petroleum to satellites now
follows this pattern of reliance on the computer. Moreover, the needs of
national defense have posed scientific and engineering design problems as
difficult as any ever encountered. Numerical computation and applied
mathematics have played an essential role in dealing with such problems.
In fact, numerical fluid dynamics was born during the 1940s for the pur-
pose of assisting in the design of nuclear weapons, combat aircraft, and
conventional ordnance and is now applied widely by industry.
Most problems of engineering or scientific interest are too complex
to be modeled and computed exactly. Instead one considers a series
of approximate models and computations, each of which illuminates a
different aspect or idealized portion of the overall problem. When used by
a skilled engineer or scientist, these mathematical models greatly enhance
the judgment that goes into design decisions and reduce the amount of
expensive laboratory and field testing required. These advantages account
for the widespread use of these models.
More specifically, mathematical models are used in engineering design
problems in the following modes:
1. To provide the first test of a new design idea. Beyond common
sense and simple hand calculations, the computer model is usually the
cheapest and fastest test of an idea. This test is applied before deciding
whether to conduct a series of experiments or to build a prototype.
2. To reduce the time and cost associated with laboratory and
field tests. Usually engineering problems contain several critical design
parameters that will have involved a certain amount of trial and error in
the search for the optimum choice. The computer is used qualitatively
(Will an increase hi parameter X improve or degrade performance - or
performance parameter Y"?) and quantitatively (Which values of design
parameters Xi,.. .,X n will optimize performance?).
3. To assist in laboratory or field tests that determine model
parameters and equations. Often the model parameters are measured
only indirectly. Thus, a computer model of the laboratory apparatus
may be needed to extract the desired information from the observed data.
Usually these models are simpler than the complex engineering models,
and their defining equations can be solved with greater precision and fewer
approximations.
4. To replace laboratory or field tests. Sometimes tests are impossible
or impracticable. For example, measurement of chemical reaction rates
at extreme conditions of pressure and temperature is very difficult, and
accumulated experience through trial and error is not adequate for solving
the problem of landing the Space Shuttle, for reasons of human safety as
well as cost.
5. To improve the education and judgment of engineers and scien-
tists using the models. The mathematical models and computer solutions
provide a vast increase in the quality and caliber of the data. Thus while
the laboratory measurement may produce some overall quantity (e.g., total
flow in and out), the computer model might yield detailed velocities and
concentrations at each point of the flow field. Because the equations are
nonlinear, it is difficult to foresee all the relevant phenomena, much less to
understand then* relative importance. Thus the model becomes for users
an experimental tool that allows them to understand a problem at a level
of detail that cannot be achieved by other means.
In summary, mathematical and computer models are used because
they are faster, cheaper, and more effective.
However, models have limitations. There are limitations in the
validity of the equations used, in the adequacy of the solution algorithm,
and in the size and speed of the computer. Also, the cost, accessibility, and
reliability of computer software and, sometimes, the cost of the computa-
tion itself can be limiting factors. These limitations in some sense define
the frontiers of science, but more specifically they define the frontiers of
applied mathematical science.
Problems of realistic interest typically involve the study of diverse
physical phenomena oA many scales of length in fully three-dimensional
settings. Though essential, experimental science in these contexts is ex-
pensive and difficult. The design of modern strategic weapons systems
epitomizes these characteristics. For example, in the design of a TRIDENT
submarine, the architecture of the vessel, the design of the missile, and
the design of the nuclear warhead all need to be modeled and integrated,
which requires many thousands of hours of computer time and stretches
available computing power and modeling techniques to their limits.
Generally, design and evaluation of new kinds of defense systems,
such as remotely controlled or robotic vehicles, will require analysis of
entirely novel systems, using parallel computer architectures among other
things. The needs of industry for technological advances are similar and
are dominated by such basic concerns as pollution, depletion of resources,
energy conservation, and efficient use of manpower. It seems safe to
say that defense and industrial needs will continue to lead numerical
computation and applied mathematics into new and challenging regimes.
Depending on the complexity of the problem and the magnitude of the
effort expended, models range from excellent to merely suggestive in their
quality and usefulness. In all cases, improvement of computer modeling is
one of the most promising avenues to improved technological performance
by OUT nation.
As one of the aims of this report, the Committee wants to show and
emphasize that in the computational approaches to most of today's press-
ing and challenging scientific and technological problems the mathematical
aspects cannot and should not be considered in isolation. There is a unity
among the various steps of the overall modeling process from the formula-
tion of the physical problem to the construction of appropriate mathemati-
cal models, the design of suitable numerical methods, their computational
implementation, and, last but not least, the validation and interpretation
of the computed results. In particular, the Committee wants to illustrate
that the steps are more often than not deeply interconnected and that
the computational process may indeed be part of the model construction.
At the same tune, there are problem areas, such as turbulence, where cur-
rent theoretical research may promise a deeper insight into an important
physical phenomenon.
In line with these aims, the report uses a "matrix approach" that
views the same problem from three different standpoints. In Chapter 2, the
traditional approach is taken of discussing a number of typical problems
leading to computational modeling from the viewpoint of the scientific
or engineering discipline in which they arise. Then in Chapter 3 certain
of these problems are touched on again, this tune from the viewpoint of
the computational and mathematical difficulties that arise in connection
with them. For example, these difficulties may involve large numbers of
degrees of freedom, different scales of time and length, or singularities of
various types. In Chapter 4 the viewpoint becomes that of the numerical
algorithms involved in the computations, such as various discretization
methods, continuation approaches, and splitting techniques.
Of necessity, many topics have been left out. The list of applications,
for example, is by no means complete, and, in fact, entire areas such
as reactor safety and reactor physics are not mentioned at all. Neither
did the Committee attempt to address all computational and mathemati-
cal difficulties nor all variations of numerical algorithms. The choice of
topics was motivated in part by energy-related considerations, the exper-
tise and interests of the Committee and its advisers, and the report's broad
purpose, whose achievement would be hindered by any attempt to be en-
cyclopedic. For the same reason, the report is not intended to be a tech-
nical summary, and this is also reflected in the fact that no attempt was
made to reference the relevant literature. The report is mainly addressed
to scientifically literate readers who know how to consult the literature
when necesary.
The Committee obtained advice and technical support from many
colleagues across the country and from abroad. A list of names of all
those who helped in this work is given in Appendix A, and the Committee
is extremely grateful for all the often extensive documentations, special
write-ups, and other comments that we received. The report is the result of
a study begun by the Committee in 1981 on computational mathematical
modeling and mathematics applied to the physical sciences with particular
reference to the needs of the Department of Energy (DOE). The prepara-
tion of the report was supported by the Applied Mathematical Sciences
Research Program of the Office of Basic Energy Sciences of DOE, and the
Committee also expresses its thanks and appreciation for this support.
As stressed in the Overview, the Committee found that improvement
of the mathematical and computer modeling of scientific problems is an
important priority for our nation. The challenge is broad, and there are
no simple remedies for current shortcomings. Accordingly, the Committee
recommends the following:
1. Increased research support for computational modeling and applied
mathematics
The technological challenges of the coming decades will impose new
tests of our abilities in computational and applied mathematics, and meet-
ing the tests will require increased research effort. As illustrated in this
report, the challenges are typically multidisciplinary in nature with applied
mathematics and modeling often in a central position.
Hence, to support research in this area, multidisciplinary teams of
adequate size to make progress on these complex problems should be
encouraged; and organizational means should be devised to facilitate their
establishment, continuity, and success.
2. Increased support for computing facilities dedicated to computa-
tional modeling and applied mathematics
Ready access to modern computing systems is essential. Lack of
equipment is the critical factor most strongly limiting academic research
in computational mathematics. There is need both for conveniently usable
local equipment and for access to large-scale computers. Local facilities
are necessary for entire problems of modest size and for such tasks as
code development, interactive debugging, test runs, and graphical analysis.
Large-scale computing is essential because of the size of many typical
problems as documented in this report. In this connection, the Committee
strongly endorses the findings and recommendations in the recent Report
of the Panel on Large-Scale Computing in Science and Engineering, Peter
D. Lax, Chairman, National Science Foundation, December 26, 1982.
S. Increased support for education and manpower development in
computational and applied mathematics
Today there are unmet manpower needs in computational and ap-
plied mathematics, as discussed, for example, in Science and Engineering
Education for the 80s and Beyond, a National Science Foundation report,
October 1980. These needs are found in industry, government laboratories,
and academic institutions. The critical challenges in this area call for a
focus on quality. Specifically, the complex interdisciplinary nature of the
problems poses special educational challenges for students and young re-
searchers, and graduate and postdoctoral fellowship support for participa-
tion in multidisciplinary teams of the type discussed above would also be
helpful.
A valuable aspect of such multidisciplinary education is the interac-
tion it creates between applied mathematicians and other applied scientists
in universities, government, and industry.
2. APPLICATIONS
2.1 HYDRODYNAMIC SYSTEMS
Hydrodynamic processes touch nearly every aspect of our lives. In a report
even many times larger than the present one, we could not possibly discuss
all hydrodynamics applications. Therefore, we have chosen a few, with
the hope that they will serve to illustrate the importance of hydrodynamic
models and to point out some of the problems that occur in developing
the models.
A large fraction of the current research in computational modeling
is done with hydrodynamics applications in mind. The great variety of
responses that we received to our requests for material for this report
attest to the diversity of these applications. Among them were applications
having to do with aircraft and wing design, both at subsonic and supersonic
speeds; global weather prediction and local phenomena such as tornadoes;
water waves and ship hull design; piping networks, such as in nuclear
reactor or power plant design; geologic phenomena, such as glacier flow
or convection hi the Earth's mantle; biological flows, such as the flow of
blood in the heart; and chemically reacting flow, such as combustion.
The general system of equations governing hydrodynamics are called
the Navier-Stokes equations. They are a statement of mass and momen-
tum conservation, the momentum equation being a formulation of
Newton's second law, F = ma. The Navier-Stokes equations were first
developed in France by Navier in the early 1800s. They represented an
improvement over the Euler equations that were first derived in 1755, in
that the Navier-Stokes equations included viscous effects that were absent
in the Euler equations. However, it was not until 1904, when Prandtl de-
veloped the boundary-layer approximations, that predictions of practical
viscous flows could be made. Practical solution of the full Navier-Stokes
equations had to await the development of modern high-speed computers
beginning in the 1940s.
Why are equations that have been known for over 100 years so hard to
solve? The answer lies largely in their inherent nonlinear characteristics.
Immediately upon looking at the equations one sees that the convective
transport terms (the acceleration in F = ma) involve velocity times its
8
gradient. This nonlinearity is always present, and it is responsible for
the existence of complex phenomena such as shock waves and turbulence.
In principle the Navier-Stokes equations alone provide a description of
turbulence; however, one would have to resolve such small length scales in
their solution that this approach is not of practical importance. Therefore,
a great many approaches to approximating turbulence effects are being
pursued. Typically, these models introduce further nonlinearities into the
system.
Other problems that arise in the solution of hydrodynamic problems
are related to the disparate time and length scales that must be resolved.
In particular, convective transport is characterized by the fluid velocity,
whereas pressure waves travel at the sound speed-typically orders of mag-
nitude faster than fluid speeds. At the same tune the effects of diffusive
processes (e.g., shearing stresses) are felt instantaneously throughout the
flow. In some cases the fluid can react chemically. In those cases, the reac-
tion rates display a strong nonlinear dependence on the fluid temperature.
This introduces yet more characteristic scales into the models.
Depending on characteristic parameters such as the Reynolds number,
the solutions to the Navier-Stokes equations either can be smooth and
steady or they can exhibit regular oscillations. Or they can be completely
chaotic. It is clear that depending on the regime the appropriate solution
procedures could be quite different.
Finally a word about boundary conditions. In some problems the flow
is enclosed, and hence boundary conditions are applied at the boundaries
of the enclosure (e.g., an automobile engine cylinder), which can often
have a complex shape. In other cases, such as an aircraft wing, the flow
is effectively unbounded, and the boundary condition should be applied at
"infinity" (see Section 3.4). Some modeling problems arise in approximat-
ing infinity by some finite boundary. In still other problems, such as
the flow of blood in a heart, the boundary is both complex in shape and
deforms depending on the forces exerted on it by the fluid. Modeling such
a problem is clearly a challenge.
2.1.1 Wings, Wind Tunnels, and Computers
The economics of the energy shortage implies that planes will fly at speeds
close to but less than the speed of sound. At such speeds there is a
"supersonic bubble" over the wing, where the local velocity of the air
relative to the whig is greater than the speed of sound. In this case the
presence of shock waves is typical and undesirable. They are undesirable
/\ /
gfS
I /
t- J
/ / /
o g
* *
l S
V
s &
** "
1 p p
J
a *
.* /
/* /
~ s?
/
3 ^
1
2 u
(
5^
\
n cd
7 1 s
1 1 ill 1 . 1 .
OH fl>
5 "^
m C^4
09' 1- MM- Of- 1 Oh'- 00-'- Oh- OS' ; l
CO t-1
a
13 "
O -4^
;| g
Q oQ
q) O
o. "^
^^
Xl \
J "CO
' Id
*
+* "
\ o3 O
}
t
a a
/
o o5
I
rt jj '*r n
/
? a e c^ 1
g > C
I :'
~. p) 2
I i
(=1 S cJ*
I i
^ : ^S^
^ 3 3
J l!S
U "rf
4? *"*
D M \
Oh 1 - 00'- Ok- U' 0-I
g qs o
dO
O m
a .a o
/
c
10
because the drag can be computed as being proportional to the third power
of the shock strength. The goal of efficient wing design is to produce
wing shapes with no shocks or only weak shocks in this transonic region.
A general mathematical theory shows that shockless wing foils exist for
given transonic cruising speeds. However, the problem of finding such wing
shapes is both overdetennined and extremely sensitive to small changes in
the data, i.e., "ill-posed" (see Section 3.6). The solution to such ill-posed
problems is still valid from an engineering point of view because operation
at neighboring off-design conditions will produce only weak shocks and
small drag, see Fig. 2.1.
Using computer algorithms created by applied mathematicians, it is
now possible to solve both the inverse problem (design) and the forward
problem (of determining the flow field for a given wing shape and velocity)
with sufficient accuracy that the use of costly wind-tunnel experiments can
be greatly reduced. This accomplishment is a striking success of recent
applied mathematics.
The theoretical areas that have contributed to this study include the
theory of nonlinear elliptic equations, complex function theory, and mixed
problems, for which a prototype is the Tricomi equation,
which is elliptic for y > and hyperbolic for y<Q. The elliptic region
corresponds to the subsonic region, and the hyperbolic region to the su-
personic "bubble." In a numerical method for the design problem, an
analytic continuation makes the elliptic region hyperbolic. The resulting
equation is solved by the method of characteristics, and then the analytic
continuation back to real values of space is performed numerically.
2.1.2 Chaos, Turbulence, and Droplets
Turbulence produces a boundary layer along the trailing edge of an aircraft
wing. The boundary layer degrades the wing performance and thus is an
important part of the design problem. The flow behind the trailing edge of
a wing contains a vortex sheet, and the roll-up of this vortex sheet produces
turbulence that constitutes a safety hazard for small aircraft flying in the
wake of a jumbo jet. The axis of the vortex roll-up is perpendicular to the
wing, and so the roll-up is intrinsically three dimensional. In the simpler
case of two-dimensional flow the vortex sheet is a line or curve in the
plane, emanating from the trailing edge of the wing foil. In regions in
11
which the line is stretching, it is geometrically stable. In regions in which
it is contracting, it is unstable. Instead of contracting, it forms a spiral
vortex structure and hence is stable.
The two-phase flow of water and steam in a cooling pipe, or of oil
and gas in an oil-reservoir production well, is also a problem in which
the geometrical instabilities of large-scale fluid motion are important.
Here an internal movable boundary separates regions of different material
properties. In some cases (a heavy fluid, e.g., water, falling into a light
fluid, e.g., air), the material interface is unstable against formation of
fingers. There is continuation of the nonlinear instabilities leading to
pinchoff of droplets and a chaotic regime (mist) that can be analyzed on
various length scales, as discussed below in the case of turbulence.
During the combustion stroke of an automobile engine, the flame
is quenched when it reaches the cold cylinder walls, and incompletely
burned fuel present in the combustion chamber at this time contributes to
pollution and to a loss of fuel efficiency. The rate of advance of the flame
from the spark plug to the cylinder walls is governed by the laminar flame
speed and the rate of turbulent mixing. Of these two effects, the second
is more important. The turbulent mixing is produced by vortices that
detach from the turbulent boundary layer at the wall during the intake
and compression strokes (see Fig. 2.2). Thus, an accurate modeling of
this problem requires an ability to treat a number of fluid singularities:
flame fronts, vortices, turbulence, boundary layers, and boundary-layer
separation.
The examples above show that singularities in fluid flow may be
geometrically unstable. When this instability occurs in a regime governed
by the scale-invariant Euler equations, the phenomenon is repeated on
all length scales and leads to chaotic solutions. Turbulence, vortex roll-
up, convection fingering, and droplet formation are examples of this
phenomenon, which we now discuss from a general point of view.
The Euler equations of fluid dynamics allow intrinsic singularities,
namely vortices, boundary and shear discontinuity layers, contact or
material interface discontinuities, and shock waves (Fig. 2.3 shows a com-
putation of the stretching of a vortex tube in a periodic inviscid flow).
Depending on the problem, special discontinuities such as flame fronts or
chemical reaction fronts (for fluid dynamics with chemistry) may occur.
Within the singularity, the Euler equations fail to be a correct description
of nature, and corrections (either parabolic effects or perhaps a more com-
plicated Euler equation with more state variables) may be required. As an
example, consider a shear layer (i.e., a jump discontinuity in the tangential
velocity component). Taking the curl of the Navier-Stokes equations, we
obtain for
12
= V X v = vorticity
the equation
+ (v V)u; (u
at
where v is the kinematic viscosity. To understand the significance of this
equation, we specialize to two dimensions. Then w is a vector in the z
direction, (w V)i>, and
is the total, or Lagrangian time derivative, so that the Navier-Stokes
equation says that vorticity moves by passive transport plus diffusion.
The extra term, (a; V)v above, induces vortex production as a result of
the stretching of vortex lines in three dimensions and is important for
considerations of geometrical stability as discussed below. In summary,
the diffusion term v&v of the Navier-Stokes equation is responsible for
the vorticity leaving a boundary or internal shear layer and diffusing into
the rest of the flow. Without viscosity there is no mechanism for vorticity
to enter (an initially irrotational) flow region. The Prandtl boundary-layer
FIGURE 2.2 Stretching of a flame by a vortical structure. Such stretching
is important for the efficient operation of engines; it enhances burning by
increasing the area of the flame. [From A. J. Chorin, Flame advection and
propagation algorithms, J. Comput. Phys. 35, 1-11 (1980),]
13
al.p te aacnoni 1 to
tip 30 "ittion 1 to 88
atap E0
"aUctlona 1 to IB
step
aic'Vlona 1 to 8
FIGURE 2.3 Successive stages in the stretching of a vortex tube in an
inviscid periodic flow. Vortex stretching is the mechanism by which energy
in a turbulent flow is transferred to ever smaller scales where it is eventually
dissipated, a, step 10, time = 0.65; b, step 20, tune = 0.88; c, step 30,
time = 1.04; d, step 40, tune = 1.21. [Prom A. J. Chorin, The evolution
of a turbulent vortex, Commun. Math. Phys. 83, 526-527 (1982).]
14
equations are a special version of the Navier-Stokes equations (scaled in the
normal direction, so that diffusion occurs only normal to the boundary).
Often the fluid singularities are geometrically unstable. They may
bifurcate in a predictable fashion, developing "rolls" (Couette flow) or
"cells" (Benard flow), or they may become irregular and highly convoluted
with a tendency toward chaos, known as turbulence. There is no scientific
reason to question the validity of the Navier-Stokes equation as a micro-
scopic description of physics even into the turbulent regime. However, its
usefulness as a description of large-scale fluid motions in the turbulent
region can be questioned, and some other description of fluid flow, such as
a random ensemble of interacting vortices, could be more effective.
There are three energy or length scales in which quite distinct charac-
teristic phenomena dominate. The smallest length scale is that on which
energy dissipation dominates. On this length scale, the Navier-Stokes
equations are the correct equations. The viscosity is large, causing velocity
fluctuations to be rapidly smoothed and solutions to be (locally) "laminar."
Since smoothness of solutions is a local property, we may conjecture that
all solutions of the Navier-Stokes equations should be smooth for all time.
This conjecture is the major unsolved problem of the energy dissipation
range. It is known that solutions with smooth data will remain smooth
for a short time, and solutions with smooth and small data will remain
smooth for all time. Both statements exclude turbulent regimes.
For weak solutions, Leray showed that the set of tunes t for which
v(x, t) fails to be smooth is a set of measure zero. Recently considerable
progress has been made in restricting the possible singularities of the
Navier-Stokes equations.
For longer length scales, viscous effects do not play a major role
and the fluid flow can be described by Euler equations. However, this
simplification gives rise to problems. The problems are not merely techni-
cal but reflect the intrinsically complex phenomenology of fluid dynamics.
The Euler equations are scale-invariant. Thus if v = v(x, t} is a solution, so
is v = v(ax, at). The inertia! range is the set of length scales dominated by
scale-invariant, universal physics. Whatever phenomena can occur (e.g.,
vortices) will be repeated on all length scales in the inertia! range. The
inertia! range is limited at the smaller end by viscous dissipation. At the
larger end, it is limited by the special boundary and initial conditions
imposed on the flow, which result in special (nonuniversal) flow behavior.
The inertial region is dominated by scale-invariant behavior. There
is a flow of energy from large-scale motions to smaller ones (the "energy
cascade"). This cascade seems plausible on physical grounds as a type of
third law of thermodynamics but does not have a rigorous mathematical
status. It can be explained as a consequence of the geometrical instability
15
of vortex lines and shear layers. As these go unstable, they generate
(smaller) new vortices and vorticity.
The energy cascade leads to a dimensional analysis of characteristic
exponents and to the Kolmogoroff 5/3 power law
for the energy distribution as a function of wave number fc. A discussion of
the experimental data in connection with the Kolmogoroff theory has the
vortices, which occur on all length scales in the inertia! range, as filling
space. Actually, it may be better to assume the contrary: vortices of
a given size fill only a small part of space. Then the smaller vortices,
which are driven by the larger ones, will occur only within the region of
these larger vortices, and in fact only within a small part of this region.
This is the notion of intermittency. It leads to the idea of a singular set
for solutions of the Euler equations, which is a Cantor set of fractional
dimension less than 3.
Intermittency leads to modifications in the Kolmogoroff exponent
and to a renonnalization group type picture of turbulence. Numerical
calculations to determine intermittency and energy cascade exponents
have been performed. The calculations start by tracking vortex lines in
a fluid flow and proceed through a sequence or renonnalization group
length scale transformations to focus on the singular Cantor set within
the solution.
In most problems, the inertia! region contains lengths too small to be
used directly in a fluid calculation. Its importance lies hi its role of fixing
parameters such as an effective or eddy viscosity, which are then used to
determine the large-scale motion of the fluid. The inertia! region is not
particularly well understood from either the theoretical or numerical point
of view.
The large-scale fluid motions are produced directly by the initial and
boundary conditions imposed on the flow. These motions are strongly
problem dependent. Numerical calculations and experiment are important
tools in their study as is the detailed analysis of simplified and idealized
flow configurations. An important theoretical question is the evolution of
initially unstable flow configurations. This question arises in connection
with the onset of turbulence and in connection with the energy cascade,
where large-scale vortices excite and drive small-scale ones.
In some problems (supercritical turbulence), the instability in an ini-
tially laminar flow is nonturbulent but arises from the bifurcation of a
fixed point. Further bifurcations lead to higher-dimensional tori, and a
general theory explains that generically the flow on the (sufficiently high-
dimensional) torus has a strange attractor as its limit set and that this
16
strange attractor is chaotic in nature. This picture has been analyzed in
the context of the Lorentz flow, which is the truncation of the Navier-
Stokes equations to include only a small number of modes. The strange
attractors found there have a Cantor-like structure. An example of super-
critical turbulence is Couette flow.
Subcritical turbulence occurs when the finite (noninfinitesimal) per-
turbation is less stable than the infinitesimal one. Then turbulence occurs
below the critical Reynolds numbers at which the linear theory shows
instability and may go directly to turbulent behavior without a discrete
sequence of nonturbulent bifurcations starting with laminar flow.
2.2 CHEMICAL SYSTEMS AND COMBUSTION
Prom the invention and manufacture of an enormous range of syn-
thetic materials (e.g., plastics) to the refining and burning of fossil fuels,
chemistry and chemical processes affect nearly every phase of our lives.
Naturally it is important to understand and control, as fully as possible,
many of these complex chemical processes. For example, we seek to find
new and better materials, to reduce costs, and to generate energy more
efficiently and with less pollution. Applied mathematics and computa-
tional modeling continue to play a valuable role in meeting these goals.
One of the oldest chemical processes harnessed by man is combustion.
The successful modeling of combustion provides an extraordinarily rich
source of challenges for the computational mathematician. Frequently,
combustion models have to incorporate all the difficulties of complicated
fluid mechanics coupled with complex chemical kinetics. The challenges
include developing algorithms to ensure accuracy and to reduce computer
tune and storage. The modeler also seeks appropriate simplifications that
take advantage of any special attributes of the governing physics in order
to gain more efficient computation. Since combustion contains a wide
range of chemical processes, we use it here as an example to illustrate
points of mathematical interest in general chemical systems.
Even within the topic of combustion there is an enormous diversity
of applications. The first topic that probably comes to mind is the model-
ing of internal combustion engines, and this is an important application.
Modeling is a part of ongoing research to design new types of engines (e.g.,
direct injected stratified charge), to improve fuel economy, to utilize alter-
nate fuels (e.g., alcohols), and to reduce pollutant formation. Similarly,
research for other combustors, such as gas turbines or power-plant boilers
17
benefit from computational models. Still there are many other important
combustion problems aside from power generation. For example, the field
of fire research is devoted to problems such as how fires spread in build-
ings and the behavior of various fire-retardant materials. An important
current topic in reactor safety is the characterization of hydrogen-air fires
such as occurred in the Three Mile Island accident. Another public safety
question deals with the problems of fire and explosion associated with
a liquid-natural-gas tanker accident. Problems of burning coal and coal
gasification are also topics of great current interest.
Perhaps the simplest chemical process from a physical viewpoint is
chemical equilibrium. At equilibrium all chemical reactions are assumed
to have gone to completion, and the species concentrations are such that
the mixture is in a minimum free-energy state. The mathematical com-
putation of the chemical equilibrium state is posed as a constrained mini-
mization problem. In combustion the equilibrium composition corresponds
to the products of combustion long after the combustion is complete.
Physically, the next step hi complication comes with the inclusion of finite-
rate chemical reactions. Here the mathematical problem is one of solving
systems of stiff ordinary differential equations (see Section 3.2 for a discus-
sion of stiffness). Models of shock tubes or flow reactors, which are used
frequently to probe fundamental questions in chemical reaction kinetics,
fall into this class of problems. The physical situation is complicated fur-
ther by the inclusion of fluid motion and heat and mass transport. In this
case the mathematical problem is one of solving systems of parabolic or
elliptic partial differential equations.
Consider the internal combustion engine as an example. What are
the things that we might hope to learn from modeling? Ultimately, we
hope to influence geometrical considerations such as combustion chamber
shape and component placement (e.g., valves, spark plugs, fuel injectors).
We also hope to provide a fundamental understanding, on the molecular
level, about how fuels are oxidized and how pollutants are formed. With
such understanding we can suggest ways to alter the combustion process
to advantage. In the past only power and size were important considera-
tions, and engine optimization could proceed mostly experimentally. Now,
however, there are too many parameters to optimize simultaneously, and
computational modeling is increasingly important.
Because of limitations hi available computational resources, two tacks
are taken in engine modeling. One is to consider mostly hydrodynamic
effects. Here the modeling of boundary shapes and component placement
is of primary importance (e.g., How should the piston face be shaped,
and where should the spark plug be placed?). Complex domains and two-
and three-dimensional effects are important. The models must incorporate
18
FIGURE 2.4 Velocity vectors computed in a direct-injection, stratified-
charge engine at a position near top dead center. The combustion takes
place in a swirling environment in a cup-like region machined into the
piston. Several vortices are seen to develop in the cup. [Prom T. D. Butler,
L. D. Cloutman, J. K. Dukowicz, and J. D. Ramshaw, Multidimensional
numerical simulation of reactive flow in internal combustion engines, Prog.
Energy Combust. 3d. 7, 293 (1981).]
turbulent hydrodynamic effects and sometimes phase-change processes,
such as fuel spray injections. The chemistry is usually simplified in these
models because it is not feasible to consider both complex chemical kinetics
and hydrodynamics on current computers. Figure 2.4 shows the velocity
vectors that result from a two-dimensional simulation of a direct injected
stratified charge engine. This is a new engine concept that is being studied
in a U.S. Department of Energy-sponsored cooperative program including
General Motors Research Laboratories, Princeton University, and three
National Laboratories.
19
In addition to the hydrodynamic aspects of engine combustion, there
are important unanswered questions about the chemistry. Therefore, the
second tack is to consider simplified hydrodynamic situations, such as
laminar flames, and treat the chemical kinetics in great detail. These
models address issues such as ignition phenomena and pollutant formation.
Figure 2.5 shows some species profiles computed in an atmospheric pres-
sure acetylene-oxygen premised flame (acetylene is an important reactant
in soot formation). This model used 30 chemical species and 103 reactions.
The results were computed using an adaptive mesh placement strategy,
and they resolve detailed structure within the flame. Note that the flame
is very "thin"- its thickness is on the order of one millimeter, while com-
bustion chamber dimensions are on the order of tens of centimeters.
An interesting aside is to note that these two approaches to combus-
tion modeling match corresponding approaches to experimental investiga-
tion. That is to say, it is usually not possible to measure or compute minor
species concentrations in complex turbulent flows, whereas it is possible
to do so in laboratory flames.
0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14
X(CM)
0.16
FIGURE 2.5 Species mole fraction profiles showing the internal structure
of an atmospheric-pressure stoichiometric acetylene-air flame. [From J. A.
Miller and R. J. Kee, Sandia National Laboratories.]
20
The inherently disparate time scales of chemical reactions, in other
words their stiffness, contribute to the numerical difficulty of solving com-
bustion models (see Section 3.2). The inherently disparate tune and length
scales for fluid transport, heat transfer, and chemical reaction are respon-
sible for the presence of steep fronts in the solutions. Also, there are many
degrees of freedom in the system of equations. The number of governing
partial differential equations is large because a transport equation must be
included for each species involved in the chemical reaction set. A system
of 30 to 50 species, involving 100 or more chemical reactions is typical even
for fuels as simple as methane. Also, for practical combustors, the model
must ultimately include complex three-dimensional geometries. Because
complete models of real combustors are too large for present computers,
an important challenge is to simplify the models (including the physical
assumptions) to a tractable level. In addition to these problems, there are
potential difficulties related to scaling. The temperatures are on the order
of 10 3 K, while some species can have important effects even when their
mass fractions are as low as 10~ 10 . Moreover, before the computation,
the peak mass fractions of the various species are usually known only to
within several orders of magnitude.
Many of the challenges of combustion modeling have been met by the
numerical-analysis community; however, many more await resolution. For
example, for systems of ordinary differential equations we understand how
to treat the stiffness that results from the complex chemical kinetics by
using stable implicit methods. However, for systems of partial differential
equations the application of these methods leaves many open questions
about how to treat the linear algebra and how to compute the error
estimates. Operator splitting methods are important in rendering the
linear algebra tractable for large problems. Stillness also occurs in low-
Mach-number flows due to very high velocity, but low-amplitude, pressure
waves. Usually, we do not care about the details of these waves, but
they can unreasonably limit the size of the time step in explicit methods.
Subcycling methods, rather than implicit methods, are sometimes used
to alleviate this problem. Subcycling is a form of operator splitting (see
Section 4.7) in which the invicid hydrodynamic equations are solved with
small time steps, while the viscous parts of the equations and the energy
and mass transport equations are solved using a much larger time step.
Another particularly important topic in combustion models concerns
the adequate resolution of localized behavior such as flame fronts. One
line of research considers adaptive meshing strategies in which a spatial
mesh network is adjusted dynamically so as to capture the local behavior
accurately. Other work considers front-tracking methods, where the flame
21
is treated as a local discontinuity and conservation equations connect both
sides of the front and predict its movement.
Combustion models must account also for fluid turbulence, an area
where even the underlying physics is not well understood. Here, computa-
tional models such as the random vortex method are proving valuable
in simulating turbulence effects, especially in the investigation of large-
scale turbulent eddies, the so-called coherent structures. We classify these
methods as "problem-dependent methods" because the physics and the
numerical model depend heavily on each other (see Section 4.9). Recently,
the random vortex method has been combined with a flame propagation
algorithm so that combustion events can be modeled. Figure 2.6 is a se-
quence of computer plots that shows the vortex velocity fields and flame
fronts as computed from a model of turbulent combustion behind a step.
2.2.1 Asymptotic Analysis
We mention here that applied analysis (in concert with computation) is
valuable for many problems hi combustion. Asymptotic methods, for
example, can actually take advantage of phenomena such as steep fronts.
They can thus be used to provide insights and to suggest approximations
that help to reduce the complexity of the numerical models and thus render
them more tractable. An important aspect of the asymptotic analysis is
the possibility of considering the dynamical stability of flames. Such work
may lead to fundamental understanding of such phenomena as the onset
of turbulence.
Asymptotic methods exploit the fact that the overall activation energy
of the chemical reaction is typically large, a well-known consequence of
which is that flame fronts are very thin. That is, chemical reaction is
only important in a thin region where the temperature first approaches
its burned value. On the unburned side of this region, chemical reaction
is negligible because the temperature is too low, and on the burned side,
it is negligible because the reaction has essentially gone to completion,
depleting the unburaed fuel. Mathematically speaking, this thin chemi-
,cally reactive region may be thought of as an internal boundary layer,
separating the unburned and burned gases. In the limit, of asymptotically
large activation energy, the boundary layer is infinitely thin, and we may
use asymptotic matching principles to connect, or match, the solutions in-
side and outside the boundary layer. The result is a flame sheet model in
which the solutions on either side of the sheet are connected by nonlinear
jump conditions that depend on local conditions at the front. Though the
22
resulting problem is still nonlinear, it represents a significant simplification
over the original problem with Arrhenius kinetics and is in fact equivalent
to the original problem in the asymptotic limit of infinitely large activation
energy.
The first studies using large activation energy asymptotics began to
appear about a dozen years ago and were noteworthy for providing for the
first time analytical flame speed formulas for a steadily propagating planar
flame. The most significant results to come out of this approach to date are
the predictions of instability and bifurcation phenomena, which describe
cellular and pulsating modes of flame propagation. These more complex
solutions branch, or bifurcate, from the basic solution (in this case, a
^7^^*"*""*r^"^^-o
T^f* - 1
FIGURE 2.6 A sequential series of computer plots displaying velocity
fields and flame fronts in turbulent combustion behind a step at inlet
Reynolds number of 10000. The combustion is stabilized by a recirculation
region behind the step. The unburned gas is a mixture of propane and air
at an equivalence ratio of 1/2. [From A. F. Ghoniem, A. J. Chorin, and
A. K. Oppenheim, Numerical modelling of turbulent flow in a combustion
tunnel, P/wl Trans, R. Soc. London. Ser. A 304, 303-325 (1982).]
23
steady, planar flame) as one or more parameters exceed critical values for
instability of the basic solution. Their significance lies partly hi the fact
that they represent intermediate modes of propagation of the flame in its
evolution from a laminar to turbulent form of propagation. Experimental
observations of these transitional modes in their pure form have recently
been reported, and the theory thus identifies critical parameter thresholds
that separate one form of flame propagation from another. A number
of bifurcation analyses using nonlinear perturbation methods have been
successful in characterizing these nonsteady and/or nonplanar flames for
parameter values near the instability threshold. However, numerical ap-
proaches may be used to great advantage in describing these transitional
solutions when parameters far exceed their critical values for instability of
the basic solution.
2.3 PLASMA PHYSICS
Controlled nuclear fusion provides a possible long-range energy source.
As currently conceived, a fusion device must contain a deuterium and
tritium plasma for a sufficiently long tune for net energy production. The
confinement may be effected either by magnetic fields or by simple iner-
tia effects. Inertia! confinement of a pellet, coupled with laser or particle
beam heating, is a distinct possibility, although the bulk of the world
fusion program centers on magnetic confinement of a plasma. The un-
derlying problems in magnetic confinement are the determination of the
equilibrium, stability, transport, and heating properties of plasmas under
realistic conditions. Numerical modeling and computation represent major
tools in this study. In recognition of their importance, the Magnetic Fusion
Energy Computing Center has been established at Lawrence Livermore
National Laboratory in Livermore, California. This Center is the third
largest scientific computing center in the United States, and its sole func-
tion is to provide computing capability to the scientists and engineers in
the U.S. fusion community. Fusion plasma physics spans a diverse collec-
tion of fields, with significant efforts occurring in the fields of Hamiltonian
particle dynamics, statistical mechanics, kinetic theories, and dissipative
and nondissipative single and multifluid models. The complexity of the
problems involved makes computation an integral part of the research and
development program for fusion.
Controlled-fusion confinement experiments have indicated that, even
in grossly stable configurations, fluctuations may play an important role in
detennining energy and particle transport. Thus, an understanding of the
24
nonlinear behavior of various plasma models is fundamental to describing
such behavior. For both ideal magnetohydrodynamics, a fluid model, and
the Vlasov-Maxwell system, a kinetic model, linearized equations have
been studied extensively. These analyses are appropriate for describing
small-amplitude deviations from a quiescent equilibrium but omit the
effects of mode-mode coupling and the onset, properties, and evolution
of turbulence. The development of a self-consistent model of plasma
equilibrium with fluctuations, stochastic particle and magnetic-field-line
behavior, and resulting transport will continue to be an important research
area.
Plasma fusion applications present many problems in which the equa-
tions of Hamiltonian dynamics appear. These systems describe single-
particle phase-space trajectories in the Vlasov-Maxwell theory of colli-
sionless plasmas. In addition, the trajectory of a magnetic field line in
a toroidal system is described by a nonlinear Hamiltonian. Thus, the
question of the existence and construction of adiabatic invariants, ex-
plicitly time-dependent or not, is of fundamental importance for these sys-
tems. Significant progress in this area might be made by combining ideas
from modern topological dynamics, numerical simulation, and perturba-
tion analysis. For example, if magnetic field lines are ergodic throughout
a volume rather than lying on closed invariant surfaces, as is given by
the Kolmorogoff-Arnold-Moser theory, then, owing to electron streaming,
the thermal conductivity within this volume will be very fast. Thus, such
questions as when ergodicity arises and what the properties of Hamiltonian
dynamics are under ergodic circumstances impact strongly on these ap-
plications.
Since computational modeling of the full three-dimensional plasma
equations is out of the question, current work utilizes a range of different
compromises, simplifying to different degrees either the physics or the
geometry to obtain a number of computationally tractable problems, each
of which illuminates in a distinct fashion, a different aspect of the full
plasma problem. In spite of these difficulties, the computational approach
is a major route to progress on the problem of controlled fusion. The
reason is simple. Experiments are expensive and must be supplemented to
the maximum extent possible by theory. The theory is highly complex and
nonlinear and is obtained by a combination of numerical experiments and
physical intuition. Clearly improved computational methods are one of the
methods through which progress is achieved in this area. Development of
more powerful computers will also be required, as it is hard to believe that
smart algorithms will by themselves suffice.
25
2.4 PARTICLE PHYSICS
The problem is to find the equations that describe the elementary particles
of subnuclear physics. This problem is among the most difficult to have
been considered seriously by science hi this century. In recent years there
have been significant gains in OUT understanding of the mathematical
structure of the equations of quantum field theory. An analysis of the
mathematical existence question for the quantum 4 equation
shows that in space-time dimension d < 3, the theory exists, whereas in
d > 5, there is no such theory (as would be given by standard methods)
except for the trivial case X = 0. The method of proof suggests nonex-
istence for the physical case d = 4 as well and offers confirmation from
the side of exact mathematical analysis of ideas advanced by theoretical
physics. This mathematical confirmation is much more compelling than
any confirmation yet offered from experimental physics.
The < 4 model is "near" to the Sine-Gordon equation, which is com-
pletely integrable. Among recent results are the discoveries of a large
class of stable, localized (in space), periodic (in time) solutions (analogs
of the Sine-Gordon breather) and of the existence of breather formation
resonances in soliton/antisoliton scattering. There are many related chal-
lenging problems in nonlinear mathematics and dynamical systems theory.
The problem with the 4 equation in d = 4 dimensions is one of several
reasons for considering in its place quantum gauge fields. A standard
approach to the quantization problem reduces it to the existence question
for a singular non-Gaussian functional integral over an mfinite-dimensional
space. In the case of gauge fields the gauge potentials form an infinite-
dimensional afGne space A on which the group Q of gauge equivalence
acts. The functional integral is only defined over the mfinite-dimensional
manifold
This manifold is not flat. The integration over M has been shown by exaxit
mathematical analysis to require introduction of coordinate patches. Only
locally, within a single coordinate patch can M AfS be regarded as an
open subset of Euclidean (Hilbert) space.
One of the methods proposed for understanding integration over M is
to understand the critical points A 6 M. This question leads to a study of
classical solutions of the Yang-Mills equation and a reduced form of this
equation, called the self-dual Yang-Mills equation. Here we are looking at
26
a specific nonlinear elliptic equation in d = 4 dimensions. A remarkable
analysis has led to a complete classification of its solutions, using methods
of geometry, topology, and algebra (fiber bundles, the index theorem, and
algebraic varieties).
The development of numerical methods for the study of quantum
fields is in its infancy, mainly because the problem is to a large extent out
of range of present-day computing machines. Among the methods used
on this problem, we mention Monte Carlo integration over a space-time
discretized version of M = AfQ. The discretized problem is called a lattice
gauge field, and even after discretization, the dimension of M is too high
to allow evaluation of integrals by direct quadrature. A second approach is
to generate series coefficients from cluster expansions of statistical physics
and to reconstruct the desired function space integral from a Fade* analysis
of the coefficients. Other attempts have been based on finite elements and
on the renormalization group. It may be that a collaboration of numerical
analysts with mathematical or theoretical physicists on this problem would
be beneficial.
2.5 CONDENSED-MATTER PHYSICS
2.5.1 Statistical Physics
In this subject, the equation of state, transport coefficients (pressure,
viscosity, and thermal conductivity, for example), and other macroscopic
properties of matter are related to and derived from intermolecular forces.
The area has diverse and important applications, ranging from metallurgy
to polymer chemistry to semiconductors and is an active area of research
from the points of view of theory, numerical methods, and experiment.
Here we limit ourselves to the two mathematical aspects of this subject:
mathematical theory and numerical methods.
The equations of statistical physics involve a large or infinite number
of degrees of freedom, and so the mathematical theory of use here is
analysis over infinite-dimensional spaces. Almost the same mathematical
structure arises as in the study of quantum fields, the relationship between
the two being that a quantum field is a continuum limit of a statistical
physics (crystal or lattice) model.
In a small class of models, including the two-dimensional Ising model,
exact solutions are known. A larger but still restricted class of models has
been analyzed mathematically with respect to qualitative behavior. One
27
3.43
FIGURE 2.7 Monte Carlo simulation of a kinetic Ising model compared
with x-ray and neutron deft-action measurements on an alloy of 60 percent
gold and 40 percent platinum that has been heated and then quenched to
60 percent of its critical temperature. The abcissa is a reduced momentum
transfer, and the ordinate is a reduced scattering intensity. (From M. H.
Kalos, New York University.)
issue recently addressed hi this work was the stability of surfaces, as is
relevant to the problem of crystal growth. Another is the effect of aperiodic
or random crystal structure on bulk material properties. A striking recent
development was the mathematical demonstration of the dipole binding
(Kosterlitz-Thouless) phase transition.
The numerical methods for statistical physics are basically of three
types. In molecular-dynamics calculations, one takes a large number of
particles and follows their motion by integration of the ordinary differential
equation
dt 2 l ' " ''
defined by the mtermolecular forces F*. This method is "exact" in its
treatment of mtermolecular forces, to the extent that they are known
and that quantum effects can be neglected, but it is approximate in its
treatment of statistics, since there are computer-dictated limits on the
number N of particles that can be included.
28
Next we mention the method of series expansions, applied to the
calculation of the equilibrium distribution dw and the partition Z, which
typically has the form
and
z =
where V is the intermolecular potential energy. Then thermodynamic
functions such as pressure emerge as derivatives or logarithmic derivatives
of Z with respect to the parameters (e.g., in V). If the potential energy
is a sum of pair potentials U then the identity
-4.50
-5.00
-5.50
-600
-6.50
-7.00
-7.50
-8.00
0.280 0.320 0.360 0.4OO 0.440 0.480 0.520
FIGURE 2.8 Energy of liquid helium as a function of density. The broken
line and points with errors show results of the Monte Carlo quadrature of
the many-body nonrelativistic Schrodinger equation. The solid line is a fit
to the experimental data. (From M. H. Kalos, New York University.)
29
substituted in the definition of du above yields an expansion that is rapidly
convergent, where
\e- u -l\l
i.e., where L/~0. This is the region of noninteraction, and so series ex-
pansions are especially useful to give weak corrections, for example, real
gas corrections to an ideal gas. The expansions appear to converge up
to critical points and with considerable work have been used to extract
information on the equation of state in that region.
Monte Carlo methods are used in the quadrature of very large-
dimensional spaces, such as the determination of du above, as well as in
the direct simulation of stochastic systems. These methods have become
an experimental tool of mathematical physics. In studying the qualita-
tive and quantitative behavior of a highly idealized model such as the
Ising model in equilibrium or very far from equilibrium, or of lattice gauge
theories, continuum and lattice models of polymers, and atomic models of
quantum liquids and solids, one can carry out numerical studies in which
the high-dimensional (i.e., many-body) character of the problem is not
distorted.
In application to the calculation of the equilibrium distribution du),
the essence of the Monte Carlo method is as follows: Starting from an ar-
bitrary point in a given ensemble, one modules a single particle position Xi
"at random" but usually so as to lower the potential energy V, occasionally
so as to raise V. After enough such elementary steps, convergence to the
distribution duj is obtained. As described here, the method is extremely
simple, and complications arise from the necessity to obtain convergence
in a reasonable tune for realistic problems.
Perhaps the most significant success has been in the microscopic
theory of classical fluids, where Monte Carlo modeling has provided the
"experimental" basis for the accurate expansion in \e~ u 1| mentioned
above. Another conspicuous success is the extraction of critical exponents
by joining ideas of the renormalization group to Monte Carlo simulation.
Figure 2.7 shows Monte Carlo simulation of a kinetic Ising model compared
with scattering measurements on a real alloy.
Figure 2.8 shows how Monte Carlo methods applied to the many-
body nonrelativistic Schrodinger equation give a quantitative account of
the energy of (real) liquid helium as a function of density.
30
2.6 GEOPHYSICAL APPLICATIONS
Out of a wide range of geophysical applications we limit ourselves here to
three modeling problems in connection with the discovery and production
of petroleum. Mathematically the discovery process is an inverse problem
(see Section 3.6) that of constructing geologic maps using seismic signals,
which in turn are generated by vibrations or explosions at the surface of
the Earth. As the signals penetrate the Earth, they cross layers of differing
density. This causes reflected signals to be returned to the surface, where
they are then recorded and analyzed. The raw data are very noisy, owing
to irrelevant near-surface density fluctuations. The noise is removed by
averaging signals from neighboring receptors or from neighboring source
locations. Then multiple reflections must be subtracted, and a compen-
sation (called normal moveout) must be introduced for effects of nonver-
tical signal propagation. Reflection from nonhorizontal layers generates
complicating shear waves as well as pressure waves. The subtraction of
multiple reflection signals can be based on Fourier analysis in a half-space
and a Wiener-Hopf factorization.
The correction for nonvertical signals can be based on a reduced
Helmholtz equation, (-A + fc 2 + v)u - 0, but in order to focus on sig-
nals moving in only one direction (either up or down), one takes first a
square root and then a power series expansion of the square root. This
process is known as the parabolic approximation and leads to the familiar
Schrodinger equation. Alternatively, the analysis can be based on ray
tracing and the solution of ordinary differential equations. The numerical
problem is to implement these steps efficiently, in view of the large amount
of data to be analyzed.
The production problem is to describe the flow of oil, water, chemicals,
and/or heat in a porous sandstone layer. The equations to be solved
are a coupled system of nonlinear parabolic equations. Generally the
equations describe mass conservation of individual species, and typically
some are nearly elliptic whereas others are nearly hyperbolic. In the
hyperbolic equations, coherent shock and rarefaction waves describing oil
banks, absorption fronts, and flame fronts may form (see Sections 3.2 and
3.3). Depending on a dimensionless number known as the mobility ratio,
heterogeneities, and geometrical effects (convergence versus divergence),
the fronts may become unstable with respect to the formation of fingers
(see Section 2.1).
Critical issues in this problem are the efficient solution of large sparse
linear systems arising from space and time discretizations and the con-
trol of numerical instabilities in the solution methods. The physical in-
31
FIGURE 2.9 Successive time steps in the movement of an oil-water in-
terface toward a producing well. The view is a vertical cross section,
with the well located on the left boundary of the computational region.
The instability is caused by three factors: by a heterogeneous conducting
channel near the bottom of the reservoir, by the converging cylindrical
flow pattern near the producing well, and by the better flow properties of
the less viscous displacing fluid (water). [From J. Glimin, B. Lindquist,
0. McBryan, and S. Yaniv, Statistical fluid dynamics II: the influence of
geometry on surface instabilities, in Frontiers in Applied Mathematics, Vol.
1, R. Ewing, ed., SIAM, Philadelphia, Pa. (to be published).]
stabilities mentioned above are controlled by the use of heat (to make
the oil flow more easily) or the use of a heavier pushing fluid (polymer-
thickened water in place of water or water and CO 2 in place of CC>2).
The third modeling problem we discuss is the process of lifting the
oil to the surface of the Earth. Because of the sudden drop in pressure,
considerable gas will typically come out of solution, and the resulting two-
phase mixture is in the droplet or mist regime mentioned in Section 2.1.
In contrast to the case of reservoir flow, even the basic equations for two-
phase flow in a pipe are not well understood.
32
2.7 METEOROLOGY
Accurate forecasting of the movement of large-scale weather patterns is
clearly an important problem. Here we mean the tracking of highs and
lows on the order of 1000 kilometers in diameter. In addition, we are
also concerned with the prediction of cloud patterns, precipitations, tem-
perature, and other elements of the weather in as much detail as possible.
Global predictions of the atmospheric flow are now made routinely at as
many as 12 levels in the troposphere and stratosphere.
By studying dynamic models of the atmosphere, estimates have been
made of the sensitivity of the atmosphere to small perturbations, often
quantified in terms of the error doubling tune. If two initial states of the
atmosphere differ by random variations of 1, then it is found that the
resultant states will differ by 2 in about 3 days and by about 4 in about
6 days. Hence, meteorologists have said that the weather is unpredictable
for periods of more than 2 to 3 weeks. But, the forecasters are still far
from being able to predict the weather accurately for even 1 week. The
best results are currently obtained by the European Centre for Medium
Range Weather Forecasting, whose 4-5 day forecasts are currently superior
to those of the U.S. National Weather Service.
The accuracy of short-range predictions is limited by four somewhat
distinct factors:
1. Initializations: (a) The accuracy and completeness of observational
data and (b) the compatibility of these data with the mathematical model,
which is, of course, a simplified representation of the atmosphere.
2. The limited number of grid points or spectral components used
in the truncated, i.e., numerical, model as a result of limited computer
power.
3. Missing or severely truncated physical processes such as cloud
dynamics, proper representation of turbulence, and radiation fields, lead-
ing to inadequacy of the mathematical model describing the atmosphere.
4. The inherent finite limits of predictability of certain types of non-
linear dynamical systems such as the atmosphere. (Compare the previous
paragraph.)
The study of problems related to the reduction of the influence of
these factors is at the heart of most research that is directly related to the
improvement of weather prediction.
The sparsity and inaccuracy of observational data [l(a)] are being
33
partially overcome by the use of weather-observing satellites. Concerning
l(b), considerable research is currently directed toward the analysis and
filtering of the observed data so as to make them compatible with numeri-
cal models. Items 2 and 3 are closely related since, in a general way, the
lack of a complete representation of all physical processes is related to the
incomplete numerical resolution of the mathematical model. Such physi-
cal processes include cloud dynamics (and the associated cloud physics),
boundary-layer dynamics (including detailed features of the terrain and
local heat, moisture, and dust sources), and detailed radioactive processes,
In item 4 a new concept, "finite limits to predictability" has emerged
from the work of E. N. Lorenz, namely the idea that the atmosphere and
perhaps certain other dynamical systems may have finite limits of predic-
tability regardless of the accuracy and detail of their initial conditions
and the accuracy with which they are computed. Lorenz concluded from
a simplified model that hydrodynamical systems with an energy spectrum
(in the customary sense) having a power-law exponent greater than 3
(e.g., -5/3) should be unpredictable in detail after a finite time. He sug-
gests a time of about 15 days for the atmosphere, whose large-scale energy
spectrum seems to have an exponent of about 2.8, with smaller scales
(i.e., less than say, 20 km) having the 5/3 power law associated with
homogeneous turbulence. This question of predictability is clearly related
to that of the stability of hydrodynamical systems and to modern theories
of bifurcation and chaos (see Section 3.5).
A large component of atmospheric research is concerned with topics
that have little direct application to forecasting. Their motivation may
lie in basic scientific inquiry or in applications to other disciplines, e.g.,
radio-wave transmission for satellite communications. They may involve
the study of fundamental hydrodynamic phenomena or other physical
processes several steps removed from specific application. In all of these
studies we find extensive use of analytical and numerical modeling of
physical systems and processes. In addition, a wide variety of statistical
approaches are used hi all areas of atmospheric research.
Let us mention a few more areas of study. In addition to the predic-
tion of large-scale weather patterns, computers are being used in limited
regions to study the development and movement of smaller-scale distur-
bances, such as hurricanes, thunderstorms, and tornadoes. Here adaptive
grid methods are used (see Section 4.2). It is hoped, for example, that a
better understanding of the details of these destructive phenomena will
lead to the possibility of altering their course and development.
On the longer time scale, where seasonal forecasts of average rainfall
and temperature are made over large regions of the globe, there are at
present several potentially useful, but unproven ideas. Here, however
34
some new mathematical and physical insight is needed in order to develop
a satisfactory "average" system of equations. This is an important open
problem.
On the still longer time scale of decades and beyond, determining
the effect of our mechanized civilization on the environment is a basic
problem. Significant progress has already been made, through the use
of simple climate models in conjunction with paleontology, astronomy,
geology, and volcanology, in efforts to understand the factors that influence
irregular alternations of ice ages and interglacial periods. These statistical-
dynamical climate models will be essential for predicting the long-term
effects of the observed increase in carbon dioxide and similar problems
that will almost certainly arise as our industrial civilization expands.
2.8 ASTROPHYSICAL APPLICATIONS
Astrophysical studies pose a wide range of problems for mathemati-
cal and computational analysis. The underlying physical theories in-
clude hydrodynamics, magnetohydrodynamics and plasma theory, radia-
tive transport theory, atomic and nuclear physics, and general relativity
theory, and hence various comments made in other applications sections
of this report apply to this area.
In recent years considerable progress has been made in theoretical and
numerical studies of a number of astrophysics! topics, as, for instance, the
study of stellar interiors, the formation of stars and galaxies, the spiral
structure of galaxies, the physics of supernovae and the evolution of su-
pemovae remnants, the formation of the solar system, and the behavior of
binary star systems, to mention just a few. The computational approaches
may range from relatively simple simulations to the numerical solution of
complex systems of partial differential equations.
As an example of the first type, simulations have been applied to
the study of star formation in a galaxy. It may be assumed that when a
massive star becomes a supernova the shock wave emanating from it can
compress the surrounding interstellar gas creating new stars. If at least
one new star is also a massive star the phenomenon can repeat, leading
to a chain reaction in the creation of stars. This is equivalent to a direct
percolation problem, and, as is typical for such problems, phase transitions
are involved. This is a nonlinear problem with a complicated structure in
space and time. Accordingly, analytic techniques are difficult, whereas
35
computational simulations are relatively straightforward and lead to the
development of realistic model galaxies.
As a second type of example we mention a problem that is amenable
to considerable mathematical analysis as well as to computational attack.
This is the question of the spiral structure of galaxies. The so-called
density wave theory approaches this question as a dynamical problem in
the form of the gravitational instability of a galactic dish with respect to
spiral modes. There are three basic approaches hi the calculation of spiral
structures of galaxies on the basis of density wave theory that may be
characterized loosely as the stellar model, particle model, and fluid model.
The stellar model represents the classical approach in the study of
stellar systems. In brief, a galaxy has a stellar component and a gaseous
component. The latter usually has a sufficiently small mass to be negligible
in first approximation. The basic equations governing the behavior of the
stellar component are the Boltzmann equation, usually in colli sionless form
when close encounters between stars are omitted, and the Poisson equation
for the gravitational potential. In a sense, these equations are much simpler
than the Vlasov equations of plasma physics, but their numerical solution
nevertheless poses challenging questions.
In the particle model, the stellar component is considered as a very
large but finite number of particles. As in our first example, the computa-
tional approach then assumes the form of a simulation process. In brief, the
motions of the individual stars are followed, and their gravitational field
is calculated by a self-consistent evaluation of the field. This technique
has provided valuable qualitative information. But the number of stars in
these models is usually of the order of 10 11 - 1 , and hence, by necessity,
numerical simulations are extremely limited in accuracy and provide only
few quantitative data for specified galaxy models.
In the third approach, the stellar system is considered to be a con-
tinuous medium. In this setting one may then study the characteristics
of wave patterns over the galactic disk and then dependence on the mass
distribution. In particular, the spiral structure of galaxies may be ex-
plained in terms of spiral wave patterns of some kind. Such stationary
wave patterns over a field of differential rotation are to be expected; in
fact, hydrodynamic waves over shear flows are known to exist for some
time especially in the form of self-excited modes. The general form of
the resulting theory of spiral galactic structures requires elaborate com-
putational techniques. At the same tune, asymptotic approaches have
'provided analytic results that have led to a better understanding of the
dynamical mechanisms despite their limitations in accuracy and in the
types of galaxies covered.
As noted earlier, these are only some examples of the many different
36
types of problems in astrophysics that relate closely with applied and com-
putational mathematics. Many of these problems involve extremely wide
ranges in length and time scales. For example, calculations of stellar evolu-
tions have to range over billions of years, while the dynamics, say, of the
supernovae phase takes place within milliseconds. Similarly, the size of a
neutron star differs by substantial orders of magnitude from the size of its
corresponding relevant gravity field. Thus, in particular, the mathematical
and computational difficulties, discussed in Section 3.2, are especially ap-
plicable here. Moreover, the wide range of the underlying physical theories
leads to many substantially different types of mathematical models, which
in turn require very different computational techniques.
2.9 STRUCTURAL MECHANICS
During the past two decades, the use of computers has transformed
large parts of solid mechanics into practical tools for a multitude of
technological developments. Sophisticated computational software is
employed throughout the nation's industries and research laboratories in
the analysis and design of structures and mechanical equipment. There
is a strong interaction between applied mathematics and solid mechanics.
Mathematical analysis has provided insight into model formulation and
the development of powerful numerical methods; and, vice versa, novel
engineering approaches have led to new research areas in applied mathe-
matics.
In the case of linear problems there exists now a relatively broad
experience in computing solutions for a range of problems concerning
the behavior of solid bodies subjected to specified loads. In general,
the computed results are considered reasonably reliable, and they have
been corroborated over a period of time by observation and practical
experimentation. There appears, however, to be a growing need for the
development of computable error estimates that can provide a realistic
check on the solution accuracy. Such a posteriori bounds have been shown
to be feasible; but their applications to large, realistic problems in solid
mechanics still requires considerable research and software development.
The general availability of economical a posteriori estimates would make
possible the consistent use of adaptive mesh-refinement techniques that
would reduce the cost of data preparation by the users and make it possible
to generate near-optimal solutions for a given amount of computational
expense.
37
The general trend in computational solid mechanics today is toward
extending the computational methodology to nonlinear problems. Sources
of nonlinearity in structural problems are (1) geometric nonlinearity due
to nonlinear strain-displacement relations, (2) material nonlinearity due
to nonlinear constitutive equations, (3) force nonlinearity due to nonlinear
stress boundary conditions, and (4) kinematic constraint nonlinearity due
to nonlinear displacement boundary conditions. The source of nonlinearity
affects the form of the resulting nonlinear equations and, hence, influences
the effectiveness of the solution techniques.
The numerical analysis of all of these nonlinear problems is not yet
at a satisfactory stage. Many computer programs for such problems exist,
but the mathematical basis for most of the methods used is insufficiently
understood, and there is little known about the accuracy of the computed
results. Moreover, in the case of nonlinear problems, few numerical com-
putations can be supplemented with sufficient experimental experience.
The situation is still best understood in the case of finite elasticity.
Even there the mathematical theory of the underlying equations is incom-
plete, and the approximation theory for these equations is generally based
on various simplifying assumptions that may or may not be valid for a
particular problem.
The state of the art in elasto-dynamic problems is in even worse
shape. It is known that multiple solutions may exist and shock waves can
develop. Moreover, not all solutions are necessarily physically relevant.
The questions of how to model such phenomena numerically and how to
determine the physically realistic solutions are as yet largely open. When
it comes to problems in finite plasticity even less is known. Although there
has been much progress in this area during recent years, no satisfactory
and complete mathematical model is available as yet. Especially, there
are profound mathematical and computational difficulties in modeling
phase changes, viscous effects, cracks and singularities, the growth of
cracks under dynamic loading, and the identification and implementation
of physically reasonable constitutive equations describing these materials.
2.10 NONDESTRUCTIVE TESTING AND TOMOGRAPHIC
RECONSTRUCTION
There is a pervasive need in technology to evaluate quantitatively the
integrity and the remaining reliable lifetime of components and structures,
e.g., from bridge girders to high-performance ceramic turbine disks. In
38
the past decade considerable progress has been made. This developing
technology is called nondestructive evaluation (NDE) to distinguish it
from older nonquantitative nondestructive testing practices. NDE presents
challenges to the applied mathematical sciences on many levels. A few of
the more important and topical problem areas are mentioned here.
In NDE applications a component is often subjected to some sort of
penetrating radiation with the aim of deducing information about its in-
ternal state from a measurement of the radiation field external to the part.
Examples include the use of ultrasonic radiation, x rays, and neutrons.
Because of its flexibility, relative cost, and safety, ultrasonic methods are
often used in NDE applications. The deduction of information about flaws
from the incident and scattered ultrasonic fields relies on the solution or
approximate solution of the inverse scattering problem for elastic waves.
In some regimes it is possible to develop and adapt imaging techniques
to the ultrasonic setting. In either case information about defects must
be obtained from band-limited noisy data. A fundamental limitation to
our current ability to utilize ultrasonic techniques broadly is our limited
understanding of the elastic inverse scattering problem in either the fre-
quency or the time domain.
In other important NDE applications, electrical currents are induced
in a material. These currents produce fields that vary depending on
whether a defect, e.g., a crack, is present in the material. The utilization
of these so-called eddy-current techniques depends on the ability to infer
information about the defects from the measured fields they produce. This
again is an inverse problem that is imperfectly understood.
Passive methods in NDS are also widely used. For example, when a
pressure vessel or aircraft is in service, a crack once formed may grow and
propagate. This will be accompanied by the release of acoustic energy.
This energy can in turn be monitored at selected sites. The problem then
becomes one of identifying and classifying the sources of these sound pat-
terns using acoustic emission studies, an increasingly important technique
in NDE. In mathematical terms the problem is the so-called inverse source
problem, which is beset with the same type of difficulties that the inverse
scattering problems possesses. Sparse, noisy data often taken at highly
nonoptimal locations are the raw information from which source charac-
teristics must be deduced.
The techniques of NDE have application in other areas, and much
can be learned in other applications fields that is valuable to NDE. For
example, there is a close connection between the need for inverse scattering
results in geophysics and NDE and acoustic imaging results in NDE and
biomedical applications. An area of great success in medical applications,
tomography, can provide useful information in selected NDE applications.
39
The greatest success of tomography has been in medical applications
such as the CAT (computer assisted tomography) scanner. Unlike ordi-
nary x-ray technique, which masks important features by superposing on
a single picture information on planes perpendicular to the beam, com-
puterized x-ray tomography provides pictures of individual thin slices
through the body. Several hundred parallel x-ray pencil beams are
projected in the plane of the slice, and the attenuation of the x rays in
each beam is measured separately and recorded. The procedure is repeated
for many different beam directions. An elaborate calculation then permits
approximate reconstruction of the x-ray attenuation density as a function
of position within the slice.
The idealized mathematical problem is the reconstruction of a func-
tion of two variables from its integrals along lines. This problem, as well
as its three-dimensional version, was solved by Radon in 1917 and later
rediscovered in various settings such as probability theory (recovering a
probability distribution from its marginal distributions) and astronomy
(determining the velocity distribution of stars from the distribution of
radial velocities in various directions). Of course, much work was needed
to adapt the Radon inversion formula to the incomplete information avail-
able in practice. Various algorithms for the numerical inversion of this
ill-posed problem have been proposed, with the present trend favoring the
so-called convolution algorithm on account of its speed and accuracy. Each
area of application has, however, its own requirements and may need a
modification of existing reconstruction algorithms or even a custom-made
one. Some algebraic methods, for instance, can easily incorporate a priori
information about the object to be reconstructed.
Recent advances in medical tomography include nuclear magnetic
resonance (NMR) tomography and positron emission tomography (PET).
In NMR strong magnetic fields are used to affect the nuclear magnetic
spin rate of hydrogen atoms. By varying the fields and their direction, the
plane integrals of the density of hydrogen can be measured and the den-
sity reconstructed by an algorithm based on the above three-dimensional
version of Radon's theorem. The technique is now regarded as competi-
tive with x-ray tomography for many purposes; and, of course, it is not
ionizing The advantage of PET over CAT is that metabolic processes
can be followed. A compound such as glucose is made using carbon-11
atoms which emit positrons. Photons resulting from the annihilation of
an emitted positron with an electron are detected by a bank of detectors
that can record coincidences. Recently, algorithms based on probabilistic
arguments have been proposed for the PET reconstruction problem
40
2.11 MATHEMATICAL MODELS IN THE BIOLOGICAL SCIENCES
The biological manifestations of the physical laws of the universe present
us with a rich variety of new phenomena that require the devlopment
of new mathematical tools and computational methods. We shall discuss
just a few examples of mathematical research in cardiovascular physiology
and neurophysiology, with the knowledge that there are many other areas
of biological sciences in which mathematics and computing are fruitfully
applied.
Blood flow in the heart obeys the incompressible Navier-Stokes equa-
tions, which, in turn, are simply a statement of Newton's laws in
differential form (see Section 2.1). The distinctive biological character
of the problem comes, however, from the moving boundaries that are in
contact with the blood. These include the muscular heart walls and the
CURVED PIVOTING DISC VALVE EXP.4: 1-9
STREAMLINES
\
FIGURE 2.10 Computer-generated plots showing the predicted opening
movement of a curved pivoting-disk valve mounted in the mitral (inflow)
position of the left ventricle. The curvature makes the valve open more
widely than a straight valve pivoted at the same point. It also helps
to prevent stagnation in the smaller opening of the valve. [Prom D. M,
McQueen and C. S. Peskin, Computer-assisted design of pivoting-disc
prosthetic mitral valves, J. Thorac. Cardiovaac, Surg. (in press),]
41
1.4
FIGURE 2.11 Transmission and "reflection" of a nerve impulse at a junc-
tion where the diameter of the neuron suddenly increases. The plots show
computed voltage as a function of tune at equally spaced positions. The
junction is at x = XQ, and the ratio of diameters there is 2.5:1. Note the
increase in propagation speed for x > XQ. A reflected wave is set up when
the larger fiber re-excites the smaller fiber after the refractory period of the
smaller fiber has elapsed. [From S. S. Goldstein and W. Rail, Changes of
action potential shape and velocity for changing core conductor geometry,
Biophysical J. 14, 731 (1974).]
elastic heart valve leaflets. The motions of these boundaries are not known
in advance; they must be computed along with the motion of the fluid.
These considerations have led to the development of a computational
model of the left heart that can be used in the computer-aided design
of artificial heart valves. In this model, the fluid equations are solved
by finite-difference methods on a regular, square mesh (see Section 3.4).
The boundaries are represented in Langrangian form as a collection of
moving points. Coupling coefficients between boundary markers and fluid
mesh points are computed with the aid of an approximation of Dirac's
^-function. This computer model (Fig. 2.10) has been used in the design
of prosthetic heart valves to remedy problems of stagnation and blood
clotting in the smaller opening of the valve. The model has also been
helpful in the study of disease processes, providing, for example, a possible
explanation for mitral valve prolapse.
42
Just as the flow of blood in the heart is ultimately governed by
Newton's laws, the conduction of electrical signals along nerves is ul-
timately governed by Maxwell's equations (Figs. 2.11 and 2.12). Here
also, nature provides a peculiar boundary condition that leads to entirely
new phenomena. In this case the boundary condition comes from certain
voltage-dependent channels located in the nerve membrane. These intro-
duce a nonlinearity, and the equations of nerve conduction take the form
of nonlinear diffusion equations the Hodgkin-Huxley equations. Without
FIGURE 2.12 Computed electrical signals at the input end of a neuron. A
brief pulse of current is applied at the periphery of a tree, and the result-
ing voltages are computed (logarithmic scale) at the input (BI), at succes-
sive branch points (P, GP, GGP), and finally at the cell body (SOMA).
[Modified from J. Rinzel and W. Rail, Transient response in a dendritic
neuron model for current injected at one branch, Biophysical J. 14, 759
(1974).]
43
the nonlinearity, signals introduced at one end of the nerve would decay
rapidly; because of the nonlinearity such signals evolve into a specific
waveform that propagates at a constant speed without distortion. This
waveform is called the nerve impulse, and it is the basic unit of long-
distance communication in the nervous system.
Currently, there is a wide range of mathematical, computational, and
physiological research activity related to the Hodgkin-Huxley equations.
Mathematically, there is extensive research on the basic theory of non-
linear diffusion equations. A particularly fruitful approach here has been
the use of piecewise linear models that expose the basic structure of the
equations. Singular perturbation methods have also been useful because
the equations exhibit a disparity of time scales.
An important physiological enterprise is the modification and applica-
tion of the Hodgkin-Huxley equations to other excitable tissues. In the
heart, for example, equations of the Hodgkin-Huxley type describe the
electrical processes that generate the cardiac rhythm and coordinate the
heartbeat. Mathematicians are just beginning to use these equations as
a basis for a theory of the abnormal rhythms of the heart, of which the
most serious is ventricular fibrilation. This theory has connections with
recent work on chaotic dynamical systems and transition to turbulence: it
appears that fibrilation is directly analogous to turbulence. This work has
enormous practical significance, since the principal cause of sudden death
following heart attacks is ventricular fibrilation.
Progress has also been made in the modeling of the input to the
neuron (whose output signal is the nerve impulse). The neuron integrates
information received through a tree of dendrites in which the signal-
ing mechanism is often described by the linear diffiision equation with
leakage. Mathematical modeling of the dendritic tree has had a substan-
tial impact on experimental neurophysiology. One reason for this is that
dendrites are too small to be penetrated with microelectrodes. Thus the
neurophysiologist can only record voltage or inject current at the cell body
and is forced to rely on the theory to indicate the significance of these
measurements with respect to activity in the dendritic tree.
Some major successes of the theory are as follows:
1. Elucidation of the dramatic differences between effects of a synapse
close to the cell body and effects of a similar synapse far out in the dendritic
tree.
2. A possible explanation of the role of dendritic spines in learning
and memory.
3. Prediction of the existence of dendro-dendritic synapses based
44
on a mathematical model of field potentials in the olfactory bulb. Such
synapses, previously unheard of, were subsequently found in electron
micrographs. This work led to a fundamental new concept of local in-
formation processing in an "analog" rather than a "digital" mode, i.e.,
without nerve impulses. Such processing is important in neural networks
such as the retina. Indeed, an extensive mathematical theory of the retina
has been developed, and this is another exciting area of current research.
2.12 ELECTRONIC COMPONENTS
The design and fabrication of modern integrated circuits is a complex
process. The number of devices that one can put on a chip depends on
the size of the chip and how small one can make its features. Over the
years the largest increase in the number of devices on a chip has resulted
from the continuing reduction in feature size and with this a reduction in
device size.
Consequently, process and design engineers have had continually to
redesign the process steps and then recalculate the resulting device charac-
teristics to ensure good electrical behavior. This has had to be computa-
tional because the equations are mathematically intractable and a trial and
error approach is prohibitively expensive and tune-consuming. Moreover,
experimental techniques tell us only what happened, not why. Effective
device design depends on determmmg both the what and the why by vary-
ing the problem parameters in the computational model.
The mathematical models on which the theory of semiconductor
devices rests are differential equations that describe the flow of current
(holes and electrons) under the influence of electric fields. When feature
sizes were large, the devices could be treated as though they consisted
of plane surfaces and edge effects could be neglected. This allowed the
development and successful use of one-dimensional analytical models ob-
tained by solving a system of three coupled, nonlinear ordinary differential
equations. As device sizes shrunk, these models became more complicated
and less accurate as edge effects became more important. In very large-
scale integration where device dimensions have reached a few micrometers,
these models are no longer adequate, and the coupled, nonlinear, partial
differential equations must now be solved in two and sometimes three
dimensions. These differential equations consist of a nonlinear Poisson
equation that describes the potential of the electric field and two non-
linear transport equations that describe the motion of the holes and the
45
electrons. Moreover, the smaller sizes have made some physical effects
important whereas they were previously neglected. This results in further
complications to the basic equations.
Partial differential equations play an increasingly important role
in simulating the fabrication processes. The transistors in a chip are
formed by implanting certain dopant ions into selected areas of the chip.
Subsequent high-temperature processes, such as growing an oxide layer,
will cause these atomic impurities to diffuse. Their final distribution is an
important factor in determining device characteristics. Insights into these
fabrication processes are especially important for increasing the yield of
reliable devices, which is a critical factor in their economic viability.
The design of the overall circuit to be placed on a chip leads to large
systems of nonlinear differential equations that need to be solved numeri-
cally. Then the efficient layout of the circuit on the chip introduces com-
binational and graph theoretical problems, which again pose formidable
computational problems.
46
3. COMPUTATIONAL AND MATHEMATICAL DIFFICULTIES
Computational mathematics and modeling is a relatively young science,
one that is expanding rapidly. The success and importance of the field
stems from the fact that its application provides the possibility of tack-
ling significantly more complex and difficult problems than would other-
wise be possible. Perhaps the greatest opportunity provided by a com-
putational approach is that it opens the wide realm of strongly non-
linear phenomena to systematic, relatively accurate, and efficient model-
ing, improving the chance that important phenomena can be isolated
and analyzed. Nonlinearities pervade nearly all aspects of applied mathe-
matics, and to a large extent these nonlinearities are responsible for the
difficulties that are encountered in computational modeling. Our purpose
in this chapter is to explore the source of some of these difficulties. In
Chapter 4, we will discuss the direction of some of the computational
research needed to resolve the problems.
3.1 DEGREES OF FREEDOM
There are several reasons for the degrees of freedom in a model, and
hence the size of computational problems, to increase. One is that we
attempt to increase the accuracy and complexity of our representations of
the physical conservation laws. Increased degrees of freedom come either
from increasing the number of dependent or independent variables. For ex-
ample, in chemical models the number of dependent variables is increased
by increasing the number of chemical species considered. An obvious
need for increased independent variables comes from the need to represent
phenomena in two and three spatial dimensions. However, even higher-
dimensional problems arise when the independent variables are not the spa-
tial coordinates but are various state descriptors; such higher-dimensional
problems are common in physics and chemistry, Unfortunately, the direct
application of numerical methods that work well in one or two dimensions
often are not usable in three dimensions. Therefore, increasing the degrees
of freedom may require significantly different algorithms.
47
Another way to increase the degrees of freedom in problems is to
allow the possibility of multivalued solutions. Certain physical systems
allow the solutions to bifurcate, a phenomenon discussed in Section 3.5.
Other systems have hysteresis or "history-dependent" properties. In these
problems the solutions depend not only on the boundary conditions but
also on the path that the transient solution follows. Structural dynamics
models for plastic materials in which the stress depends not only on the
strain but on the straining rate as well exhibit this behavior.
Consider an example to quantify the magnitude of the computational
requirements for multidimensional problems. First suppose we want to fol-
low the evolution of a chemical model having N chemical species. Suppose
also that the problem is stiff (see Section 3.2) so for each time step in
the evolving system we solve N nonlinear algebraic equations in N un-
knowns. .Even for models with a large number of species this is a relatively
straightforward task. However, if we now want to introduce transport
phenomena, such as fluid mixing, the problem has to include the spatial
dependence of each chemical species.
For a one-dimensional case assume that we use I spatial mesh points
and that we estimate what each species will do on the basis of its current
local value and that of its immediate neighbors (a three-point spatial stencil
involving all species). We now have NI unknowns and a nonlinear system
that has 3NI nonzero entries in each equation. It is typical to have 40 to 50
species and to require 100 mesh points to resolve the species concentrations
accurately. In this case over three quarters of a million words of memory
are needed just to store the approximating local linear system along with
the solution. This alone is larger than the fast-access memory of most
modern computers. Suppose further that the same model is to be posed
in two dimensions on a I x I mesh and in three dimensions on an / X / X /
mesh. The solutions themselves require NI 2 and NI 3 words of storage,
respectively. Worse, the totality of coefficients in the typical linear system
will be 5N 2 / 2 and 7N 2 I 3 . It is of practical interest to want solutions to
three-dimensional combined kinetics and transport problems. However,
even for a modest system of 20 species and 50 mesh points per spatial
dimension, 2.5 million words are needed to store the solution alone. An
additional 350 million are required to store the coefficients of the linear
system. (For comparison, we note that the largest computers currently
available have 4 million words of fast-access memory.) Moreover, in the
high-dimensional cases, the linear system is not conveniently structured
for efficient solution.
Clearly a major problem in solving high-dimensional systems of par-
tial differential equations is that after discretization the resulting system
of approximating linear equations can be much too large to be solved
48
effectively by direct means. Except in simple cases the resulting equations
have to be solved iteratively, and, especially for strongly nonlinear equa-
tions, the convergence properties of the iterative process may be a major
concern.
Operator splitting methods, such as alternating direction implicit
(ADI) (see Section 4.7), also provide a current effective approach for
high-dimensional problems. The approach here is to alternate the solu-
tion of a set of lower-dimensional problems. The alternated problem ap-
proximates the original problem with sufficient accuracy, but the set of
lower-dimensional problems is much more easily solved, even in aggregate.
Again convergence and accuracy properties must be established in all but
the simplest of cases.
For some problems it may be more efficient to depart from the con-
ventional ideas of discretization on a mesh network and consider instead
mesh-free methods or (see Sections 4.8 and 4.9) Monte Carlo methods.
Unknowns can be changed: instead of the amount of material at a given
location, one can ask for the amount of a given wavelength in the solution
as a whole. In fluid mechanics, vortex methods could be a more efficient
approach. Such methods reduce the size of the linear algebra problems, in
comparison with the mesh-oriented methods, if the new unknowns carry
all the important information desired.
Decisions about which methods are most effective may depend
strongly on advances and changes in computer architecture. Since new
architectures, such as vector and parallel processors, are now evolving,
the numerical analyst has to re-evaluate his approaches periodically. Also,
considerations such as the relative cost of memory versus central processor
time can bear heavily on decisions regarding algorithms.
3.2 DIFFERENT TIME AND LENGTH SCALES
In principle a model can contain important length scales that range from
the size of an atom to the size of the universe. In practice, however, we
limit the range of scales by approximation (see, for example, Section 3.7).
Nevertheless, it is often the case that mathematical models of physical
processes are characterized by the simultaneous presence in their solution
of significantly different time and length scales. The solutions to these
models will have regions of strongly localized behavior, such as shocks,
steep fronts, or other near discontinuities. Therefore important topics of
research in numerical analysis are the consideration of such circumstances
49
and the development of efficient theories and methodologies for their com-
putational solution. Indeed, we often find situations for which solutions
are not possible, or at least not practical, without the application of spe-
cialized methods to deal with the multiple characteristic scales.
Fluid mechanics and chemical kinetics are two areas that provide a
rich source of examples for multiple and disparate scales. Fluid-mechanical
processes are commonly characterized by groups of parameters that are
indicative of the various scales in a problem. Some examples are the Mach
number (relates velocities to sound speed), the Reynolds number (relates
inertia! forces to viscous forces), the Prandtl number (relates viscous effects
to thermal effects), and the Damkb'hler number (relates chemical reaction
rates to diffusion rates). When any of these numbers is very large (or
small) it is likely that the solutions to the models will have regions of
localized behavior. For the case of large Mach number the possibility of
shocks exists. Similarly for large Reynolds number we expect to encounter
boundary layers in the vicinity of solid walls. When the Prandtl number
is large we expect that thermal boundary layers will be much thinner
than viscous boundary layers. And when Damkohler numbers are large we
expect narrow reaction fronts.
In chemical kinetics we find large numbers of chemical reactions com-
peting for the available chemical species at widely different rates. As a
result some species are either consumed or produced rapidly or slowly,
while other species are being both produced and consumed simultaneously
at high rates, with their net production rate being relatively slow. This
chemical behavior is responsible for the many widely differing time scales
in the mathematical models. The computational models of these processes
are characterized as either multirate problems or stiff problems.
It is worthwhile to point out the distinction between multirate prob-
lems and stiff problems. In both problems the system itself is equally
capable of rapid or slow changes. Multirate problems are those in which
at least one component of the solution is changing rapidly, even though
others are changing slowly. Numerical methods for these problems must
take time steps that are small enough to resolve the fast transients, so they
are controlled by accuracy not stability considerations. Stiff problems, on
the other hand, are those in which all components of the solution axe
changing slowly compared with the fastest characteristic scales possible in
the model. In these cases explicit numerical methods are forced to take
much smaller tune steps than are needed to maintain accuracy in order
tomamtam stability. Often problems that begin as multirate problems
become stiff problems as an equilibrium or steady-state condition is ap-
proached. Stiff problems are usually solved efficiently by implicit methods.
50
For stiff or multirate problems, it is perhaps useful to consider methods
that treat the fast components differently than the slow components.
Given that multiple-scale problems are of practical importance, we
must consider: why are there computational difficulties, and what can be
done to ameliorate the difficulties? In the case of disparate length scales
one difficulty is that of representing accurately the highly localized be-
havior. If the solution is represented discretely on a mesh network, the
mesh must be sufficiently fine to capture the localized behavior accurately.
The whole topic of adaptive meshing is critically important for these prob-
lems (see Section 4.2). Here, instead of computing on a fixed prespecified
mesh, the mesh adjusts itself dynamically as the solution develops hi order
to maintain accuracy in the solution.
For situations in which the localized behavior is known to be ap-
proximated well by very sharp fronts (e.g., shocks or flames), front track-
ing methods can have significant advantages. Unlike the adaptive mesh
approach where the solution is resolved smoothly through the front, the
front tracking methods approximate the front by a discontinuity whose
magnitude, speed, and location are to be found. Then elsewhere in the
region the conventional discrete representations are adequate.
In problems like chemical kinetics, the disparate time scales cause
the governing differential equations to be stiff. Here, explicit solution
methods are well known to be extremely inefficient, and some form of
implicit method is needed. For systems of ordinary differential equations
the problem has been worked out, and high-quality computer software is
available. However, when stiffness is encountered in the context of systems
of partial differential equations the remedies are much less developed. The
same techniques used for ordinary differential equations, when applied
directly the partial differential equation problems of practical interest,
often yield problems that are simply too large for current computers.
Several approaches show promise. One is to develop operator- splitting
methods in which the stiff parts of the problem are split off and solved as
a series of smaller and hence more tractable problems. Another approach
is to attempt to remove the stiffness by solving instead an approximate
(yet sufficiently accurate) system of equations. This tack benefits from an
asymptotic analysis of the equations. Both approaches have found recent
successes in fluid mechanics and in combustion chemistry.
51
3.3 SINGULARITIES IN COEFFICIENTS, DATA, OR STATES
Difficulties similar to those encountered in multiple-scale models are often
found in problems having singularities in coefficients or states. That is, a
singular or discontinuous coefficient can give rise to localized behavior in
the solution, such as very steep fronts. An example could be a material
interface in a structural or heat-transfer problem, say between a steel and a
plastic part. At this material interface the solution (stress or temperature
gradient) might change rapidly. In order to maintain accuracy in the
computed solution, the numerical procedure would have to resolve this
frontlike behavior. The situation is analogous to the occurrence of a shock
or a flame front. However, we usually know where the material interfaces
are, so they are perhaps easier to deal with than phenomena such as shocks.
Phase transitions can also produce discontinuous coefficients. Take
a model in which a melting front is traveling through a region. Usually
the properties of the molten material are quite different than those of the
solid material. In fact, different governing equations may even be required
for the two regions (e.g., fluid motion may be modeled in the liquid but
not the solid). In any case, the solution is likely to exhibit jumps in its
properties at the transition, and the numerical method will need to locate
and resolve it. This situation is more like a shock, in that the position of
the phase transition front is not known a priori, and thus the numerical
method must both locate and resolve it.
3.4 BOUNDARY CONDITIONS
The solution to a boundary or initial-boundary-value problem depends
strongly on the boundary conditions. Thus it is important to understand
the relationship of the boundary conditions to the differential equations
and to their discrete representation. Most important, the boundary con-
ditions must be chosen so that the problem is well posed. For a large class
of problems there is a satisfactory theory of admissible boundary condi-
tions, but for many problems, those involving coupled hyperbolic-elliptic
systems or disparate time scales, for example, only a rudimentary theory
is available.
A possible error in prescribing boundary conditions for hyperbolic
equations is to overspecify or underspecify the number of boundary condi-
tions. Overspecification usually causes nonsmooth solutions with mesh os-
52
dilations near the boundary. Underspecification does not ensure a unique
solution, and the numerical solution may tend to wander in steady-state
calculations. In either case, the results are not accurate and one should
be skeptical of even the qualitative behavior of the solution. It should
be noted that the way in which boundary conditions are specified for the
difference equations can change a well-posed continuous problem into an
ill-posed (unstable) discrete problem.
Two of the most common methods used to incorporate boundary
conditions into discrete equations are the extrapolation and uncentered
differences methods. In the extrapolation method, the domain of the prob-
lem is extended and the solution is extrapolated to fictitious points outside
the integration region. The nonphysical solution at these points is defined
so that the discrete equations are consistent with as many relationships
as can be derived from the boundary conditions and differential equa-
tions. The extrapolation formula can do this best by incorporating the
discrete boundary conditions into the extrapolant. Additional relations
can be generated by differentiating the boundary conditions with respect to
time, replacing all time derivatives by space derivatives using the original
differential equation, and discretizing the resulting equations.
The uncentered differences approach is to extend the number of
boundary conditions so that all components of the solution are defined
at the boundary. Again, these additional boundary conditions must be
consistent with the original problem and as many relationships as can be
derived from it. An uncentered difference approximation is then used to
approximate the spatial derivatives at the mesh points nearest the bound-
ary.
Irregular domains can be imbedded in an underlying regular grid
that is not aligned with the boundary, or an irregular grid can be con-
structed that conforms to the boundary. The discrete approximations to
the equations away from the boundary are much simpler on the regular
imbedded grids, but the boundary conditions are difficult to approximate.
Boundary-fitted grids can be generated algebraically in the original physi-
cal domain, or the domain (and hence the grid) can be mapped onto a
regular grid in a simpler domain and the equations solved there. The
algebraic-grid-generation methods have the advantage that the equations
and boundary conditions are unchanged, but the differential operators are
more difficult to approximate on the nonuniform grid. When using the
mapping method, the differential operations are easily approximated but
the transformation can greatly complicate the equation and sometimes
obscure important properties such as the conservation laws expressed by
the equations.
53
3.4.1 Boundary Conditions at Infinity
Many physical problems require the solution of partial differential equar
tions on some infinitely large domain fL For computational reasons this
domain is often replaced by a finite domain QI. Then the difficult problem
of specifying boundary conditions at its finite artificial boundary B arises.
It is especially important that these artificial boundary conditions do not
introduce spurious phenomena. Consider, for example, a nonviscous fluid
that at subsonic speed leaves HI through the boundary B. There is one
characteristic direction that points back into the region HI, and there-
fore one boundary condition has to be specified on B. But, in general,
no detailed knowledge of the solution on B is known and therefore other
principles have to be applied. For example, if one has solved the problem
by difference approximation then one predicts the solution on B from in-
side wi, for all the dependent variables. Thus this procedure amounts to
overspecification of the solution on B. Another principle has been proposed,
namely, to specify the boundary conditions on B so that no reflection of
high frequency takes place. However, numerical experiments have shown
that such approaches do not always work.
3.5 BIFURCATIONS AND CHAOS
3.5.1 Buckling and Collapse Behavior, Bifurcations
In general, the equilibrium equations of a mechanical structure involve a
finite number of parameters, that is, they have the generic form F(x, p) =
0, where x varies in some state space X and p 6 R m represents a parameter
vector. Thus, in general, the solution set {(x,p) 6 X xR m ; F(x,p) = 0} is
a manifold X x/?" 1 , and one topic of interest is the location and character
of the singular points of the solution set.
Without going into detail, suppose that we follow some curve on
the solution manifold defined by some combination of parameter values
with one degree of freedom represented by a single variable X. Then
we may encounter certain critical points on the solution curve where
the mechanical structure may suffer a loss of stability. Such a loss of
stability actually corresponds to a dynamic phenomenon whereby the
structure undergoes a sudden change of deformation. The dynamics of this
54
phenomenon are not described by the equations of equilibrium F(x, p) = 0,
but it is possible to deduce from the shape of the (static) solution manifold
at that point what type of sudden changes may be expected.
We use a few figures to illustrate the situation. In the upper left-hand
part of Fig. 3.1 the point denoted by (1) is a so-called limit point or turning
point and an increase of the load-intensity X beyond the critical value X c
results in a jump from (1) to (presumably) (2). This type of behavior is
called a snap-through or collapse. In the case of the upper right-hand part
of Fig. 3.1 the instability phenomenon is related to the bifurcation of the
solutions at (1). This behavior is called buckling. This part of the figure
is a classical example of stable buckling where a distinct change in the
character of the load deformation is encountered when the load-intensity
X passes the buckling load X c but where no failure of the structure occurs.
On the other hand, the lower part of Fig. 3.1 shows an unstable buckling
point where we again observe a dynamic departure from the bifurcation
point to some other state [presumably the equilibrium state (2)]. The
geometrical shape of the bifurcation branch n is the determining factor in
the question of whether the bifurcation point is stable or not.
These examples already indicate that for a deeper understanding of
the behavior of a mechanical system it is necessary to analyze the shape of
the full solution manifold. Of course, the choice of the parameters entering
into the definition of the equation is of critical importance here. In essence,
catastrophe theory provides some information about the selection of par-
ticular minimal numbers of parameters, but in practice the parameters
are simply chosen to correspond to the certain natural features of the
structural problems.
In view of these observations the aim is to develop procedures for
a computational analysis of the form of the equilibrium surface. Some
methods for this purpose are mentioned in Section 4.4 on continuation
methods, but the entire problem is still a wide open research problem. In
particular, it has to be noted that we can compute only points that are
approximately on the solution manifold of a discretization of the original
problem. Thus, the questions arise whether phenomena, such as limit
points, or stable or unstable bifurcation points, encountered on the solution
set of the discretized problem actually correspond to similar phenomena
for the original problem, and, if so, what errors have been encountered.
These questions are as yet largely unanswered and represent considerable
research challenges. In this connection, it might be mentioned that for
nonlinear problems of this type the solution manifold of the discretized
equations often has a different number of connected components than
that of the original equations. The components that do not approximate
the exact solution manifold have been called spurious, or numerically
55
F=XF
dliplacamant W
cylindrical shell - roof segment
Collapse Behavior
F=XFo
dliplooamant W
simply supported plate
Stable Buckling
F=XF
dliplooamant W
cylindrical shell segment
In compression
Unstable Buckling
FIGURE 3.1 Instability behavior of shells under different loadings F. On
the right the load intensity X is plotted versus a characteristic deflection
W. [Prom W. Rheinboldt and E. Riks, "A survey of solid techniques
for nonlinear finite element equations," State of the Art Surveys of Finite
Element Technology, Chap. 16, special AMD pub!., American Society of
Mechanical Engineers, New York (1983).]
56
irrelevant solutions. They are being observed more and more often in
various applications, but their study has only recently begun.
3.5.2 Chaos in Deterministic Systems
It has been remarked that one of the most engaging problems in nonlinear
dynamics is that of understanding how simple deterministic equations can
yield apparently random solutions. As is now widely recognized, systems
with this behavior appear hi many of the branches of scientific endeavor.
This near universality has come to the attention of physical scientists and
applied mathematicians principally through the use of computers in the
study of properties of dynamical systems.
The onset of chaos in deterministic systems (the stochastic instability)
signals an unusual sensitivity to initial conditions, i.e., the trajectories hi
phase space diverge widely as tune goes on even though the initial con-
ditions are arbitrarily close to one another. The behavior is such as one
would expect hi a space with everywhere negative Gaussian curvature. On
the other hand the trajectories may hi fact tend to a single orbit nearly
filling a subspace of the phase space (a strange attractor) in an ergodic
manner. In even the simplest systems of this type it is possible for several
such attractors to coexist side by side, with the initial conditions deter-
mining which one is reached asymptotically hi tune. These and other
numerous unusual properties of dynamical systems displaying chaotic be-
havior have led to some open problems in computational mathematics.
Examples of such problems include the following: When on varying a
parameter the stochastic instability sets in, the continued use of an algo-
rithm describing the evolution of the system prior to the onset of chaotic
behavior may no longer be appropriate. In this event it may be more prac-
tical to take advantage of the stochastic nature of the system and use more
or less conventional statistical methods. In order to effect such a change in
computational methods it is necessary to detect the change from orderly
behavior of the system to chaotic behavior. Theoretically, the sensitivity
of the system to initial conditions accompanies the occurrence of positive
Lyapounov exponents, i.e., numerical indices of the asymptotic rate of
divergence of initially nearby system trajectories. Since the Lyapounov
exponents are defined as time asymptotic quantities, straightforward com-
putation of these quantities requires a considerable computational effort,
including a calculation of the evolution of two initially nearby trajectories,
or the simultaneous integration of the associated variational equations.
An important tool for studying the properties of dynamical systems
57
in the chaotic regime is the first return map, sometimes referred to as
the Poincar6 map. Such a map displays in a graphical manner various im-
portant aspects of the attractor. Unfortunately, the simplicity of generat-
ing such a map is restricted to systems no larger than three dimensional,
because of the self-evident difficulty of making multidimensional graphic
displays. As the study of higher-dimensional dynamical systems advances
there -will be an ever more urgent need for a higher-dimensional equivalent
of the first return map.
As mentioned earlier, in the chaotic regime it may be more appropriate
to describe the properties of the system trajectory in statistical terms than
in terms of a trajectory evolution. In order to obtain such a statistical
description it is necessary to have appropriate information about the in-
variant measure associated with the given dynamical system. With few
exceptions such measures cannot be derived a priori but must be obtained
from detailed calculations of the trajectories. A finite machine computa-
tion of the system trajectory will inevitably introduce some errors that
may be all the more serious in the chaotic regime because of the previously
mentioned exponential divergence of nearby trajectories. There is then the
question of how accurate the computation of a chaotic trajectory must
be in order to yield enough information for constructing the appropriate
invariant measure.
While the above examples have been drawn from dynamical sys-
tems representable by ordinary differential equations there are other sys-
tems, e.g., those described by partial differential equations and by integral
equations, whose study is likely to be replete with similar computational
difficulties.
3.5.3 Symmetry Breaking
Bifurcations commonly arise in connection with a loss, or breaking of
symmetry. In such cases the extra structure of the problem symmetry
may simplify the analysis. We give an example from fluid mechanics, A
classic experiment concerns flow between two rotating vertical cylinders.
For small angular velocities of the cylinders the flow is laminar and is called
Couette flow. As the angular velocity is increased the vertical translation
symmetry is broken and axisymmetric Taylor cells appear. These cells
resemble a row of adjacent smoke rings, each one rotating in a sense
opposite to its immediate neighbor. A further increase in velocity breaks
the axial symmetry, and the Taylor cells become wavy with a time periodic
shape. Eventually the flow becomes turbulent. Specific experiments on
58
these systems are a classic topic in fluid mechanics. However as often
happens in science, one cannot see important phenomena in the absence of
a carefully thought out theory. Recently a number of precise experiments
have been conducted concerning the role of strange attractors. Routes to
turbulence and, hi particular, wavy Taylor cells were carefully observed.
The results confirmed some of the theoretical predictions but not all of
them. One theoretical proposal was that generic limiting sets, known as
strange attractors, would dominate turbulent flow patterns. This proposal
was supported by the fact that strange attractors were observed in low-
mode approximations to fluid flows. More refined calculations with more
modes included indicate an absence of these strange attractors, but the
phenomenon is still indicative of the highly complex solution manifolds
that can arise hi nonlinear problems.
3.6 ILL-POSED INVERSE PROBLEMS
The notion of a well-posed problem is due to Hadamard: a solution must
exist, be unique, and depend continuously on the data. The term "data"
can have a variety of meanings; in a differential equation it could include
any or all of the following: boundary values, initial values, forcing term,
and even the coefficients in the differential equation. Since data cannot be
known or measured with arbitrary precision, it was felt for a long time that
real physical phenomena had to be modeled by mathematically well-posed
problems. This attitude has changed considerably in recent years, and it
is now recognized that many applied problems are ill-posed, particularly
when they require numerical answers in the presence of contamination of
the data.
El-posed problems often arise in the inversion of well-posed problems.
Consider, for instance, a well-posed problem that smooths the data or
attenuates its high frequencies. The inverse problem, in which the role of
data and solution are interchanged, will then be ill-posed. A simple but
important example is the Fredholm integral equation of the first kind
/ k(x, y)u(y)dy = v(x), < x < 1
Jo
or, in operator form
Ku = v
Assuming the kernel to be continuous on the closed unit square, the
59
Riemann-Lebesgue lemma gives
yi
lim / k(x t 2/)sin nydy =
n ^oo JO
so that K attenuates high frequencies and thus transforms widely different
functions u into approximately the same v. If v is regarded as the data, the
solution u of the integral equation, therefore, cannot depend continuously
on the data. Indeed, it is not necessarily uniquely solvable, either take,
for example,
N
If, unaware of these difficulties, one attempted to solve the integral equa-
tion by discretization, one would find that the corresponding matrix prob-
lem is singular or ill-conditioned (and the finer the discretization, the more
ill-conditioned the matrix problem).
More general inversion problems can also be reduced to an equation
of the form Ku = v, where K is a continuous transformation, v is the
data, and u is the solution being sought. It may happen that u itself is
not the quantity of principal interest but rather some functionals of u such
as some of its moments or its values at a few specified points.
Thus, there are three pieces of information that are central to the
numerical resolution of an inverse problem: (1) the model M, representing
the equation involving a mapping K between appropriate spaces; (2) the
observation operator 0, representing the measurements that can be made;
for instance, we might have O(w) = {v(xi)-,i = l,2,.,.,n}; (3) the intel-
ligence operator J, which specifies the information we wish to extract from
the solution. For instance, we might have J(u) = {fx k u(x)dx; k = 0, 1, 2}.
Most of the existing approaches focus almost exclusively on the model
M, taking account of a priori information about the solution such as
smoothness, positivity, and bandwidth. The Tychonov regularization
method is of this kind. By restricting K to a suitable subspace it may
be possible for the restriction of AT to be one to one with a continuous in-
verse. It might be more appropriate to study the triplet M,0, J for different
choices of O and J. One goal of such a new approach would be to classify
and quantify ill-posedness by developing comparison principles and order
structure.
Hi-posed problems with their attendant numerical difficulties abound
in practice. For instance, in scattering from obstacles one may wish to
determine the shape of an object (or its surface impedance, if its shape is
60
known) from far-field measurements. Coefficient identification problems
arise in many contexts, one such problem being the determination of the
sound speed in subsurface media from measurements of the field at or
near the surface. Another class of ill-posed problems is associated with
image reconstruction (tomography), the medical applications of which are
now widely used (CAT scanners, for example) (see Section 2.12). More
recently, similar ideas have been applied to the nondestructive testing of
mechanical structures to detect cracks in fuel rods, weaknesses in rotor
blades of jet engines, and faults in screws in ships. Optimum filtering
and inverse problems in Fourier optics (restoring data designed by a band-
limited filter, for example) are other areas of current research.
3.7 EFFECTIVE MEDIA
The study of "bulk" or "effective" parameters for composite media is of
fundamental importance. Depending on the particular application area it
may be relevant to consider a periodic (or an almost periodic) or a random
formulation.
3.7.1 Homogenization, the Deterministic Approach
Many problems of physical interest involve several length scales. As an
important example we mention the study of composite materials in struc-
tural mechanics. Owing to a particular manufacturing process, a distinct
structure, e.g., periodically or almost periodicity, is often present.
Homogenization is an approach for deriving the macroscopic
properties of the material from the known microscopic ones. A variation
of this is to replace a complicated geometrical configuration by a simpler
one, e.g., to replace a framework by a plate or to smooth out a rough
surface.
Homogenization may be applied to linear and nonlinear problems and
can provide qualitative information about physical macrolaws as well as
good approximations to the various parameters appearing in these laws.
The most important mathematical tool for the deterministic approach
is some form of asymptotic expansion. The theoretical results concern the
limiting behavior when certain of the relevant scales become very small.
61
One is likely to obtain different homogenized formulations from different
limiting relations between the length scales.
The numerical solution of a problem introduces a new length scale:
the mesh size in a finite-element method, the step size in a finite-difference
method, or the wavelength in a spectral method. Different limiting rela-
tions between this new length scale and the original physical length scales
usually lead to different results. In this sense the numerical treatment has
to be considered simultaneously with the homogenization process. One has
to design algorithms that will adaptively select the correct homogenized
formulation and discretize it appropriately.
Whereas much theoretical work has been done within the last ten
years analyzing the effects of multiple scales in continuous media, the
study of numerical discretizations of such problems is still very much in
its infancy.
3.7.2 Analysis of Random Media
Assessing the effect of random fluctuations in the coefficients of a par-
tial differential equation is a basic mathematical problem that arises con-
stantly in science, applied science, and engineering. Finding effective con-
ductivities of composite conducing materials such as soil or a metal alloy,
finding effective fluid equations for flow in porous media (Darcy's law),
determining the rate of sedimentation of particles in a fluid, and many
other problems are problems that may require analysis of random media.
In addition to the quantitative aspects of the problem, many interesting
qualitative questions can be posed as well. For example, what is the nature
of the spectrum of the Schrodinger equation with a random potential? It
has been shown that in one dimension, for a large class of random poten-
tials, the spectrum is always discrete. For randomness of large magnitude,
it was recently shown that there is no diffusion, without restriction on the
dimension of the space.
From the applications' viewpoint, one can frequently model ade-
quately a random medium, such as a suspension, for example by a con-
tinuum with suitable constitutive properties. The continuum equations
may be linear or nonlinear, and the constitutive laws may be known
only qualitatively from experimental data. In such a context, theoreti-
cal investigations are useful and necessary in order to understand how
the phenomenological continuum equations arise from the known micro-
scopic Jfrueta* One can then find mathematical characterizations for
the relevant constitutive laws that can lead to interesting conclusions.
62
Examples of this arise frequently as, for instance, the determination of
bounds for the effective conductivity of composites.
The effective media ideas of Maxwell in the last century dominated
theoretical calculations for a long time. In the last two decades, press-
ing technological needs have caused a major expansion in materials re-
search. In attempting to achieve a mathematical understanding of the
methodological basis for effective media calculations that abound, one finds
a lack of theoretical foundations.
Mathematical methods in random media are drawn from analysis,
probability theory, asymptotic methods, and differential equations. The
goal is to develop tools that bridge the gap between microscopic and
macroscopic descriptions, give qualitative information about constitutive
laws, and determine when residual stochastic effects remain and how they
can be characterized. In addition, it is of interest to find important specific
problems on which more detailed analysis, including numerical analysis,
can be carried out so that these problems can serve as benchmarks in the
theoretical development.
3.8 VALIDATION, ERROR ASSESSMENTS, AND
SENSITIVITY ANALYSIS
The results look good, but are they really right? The question is often
not answered satisfactorily. Indeed, often too little attention is paid to
the difficult topic of model and code validation. Once the results look
plausible, we are often either unable or unprepared to take the valida-
tion process further. The validity of a model depends both on having a
proper physical model (Do the governing equations adaquately represent
the physics?) and on having an accurate computational representation of
the governing equations. The mathematician has a responsibility in both
areas. He should help determine the "well-posedness" of the models and
from a mathematical point of view help the physical scientist determine
the appropriateness of the model. Given that the equations are proper, the
computational mathematician must be sure that the numerical procedures
used accurately approximate the solutions of the governing equations.
Frequently modeling is done in conjunction with experiments, and
those results are used to validate the model. Comparison with experiments
simultaneously tests both the validity of the governing equations in the
model and of their numerical solution. In this case the source of any
discrepancy is not easily isolated. A reasonable goal should be to validate
the numerical procedures and the physical models separately so that a
63
model can be developed more confidently. The design of reliable numerical
error estimates for the computational methods is an important step.
3.8.1 A Posteriori Error Estimates
An important current research topic is the development and application
of a posteriori error estimates. Such procedures are valuable for models in
which differential equations are approximated by their discrete analogs on
a mesh network. When a posteriori estimates of the error associated with
the discretization of the continuous model are possible and they often
are then potentially it is possible to control the errors. Adaptively moving
the mesh network to control or mimimize the error is one application of
this approach. With a posteriori error estimates, it should be possible,
in principle, to give a strict bound for the error on completion of the
computation.
The selection of the specific accuracy requirements depends strongly
on the goal of the computation. Often it is desired to obtain detailed
information about the solution itself; in other cases, the main focus is the
value of a specified functional of the solution, as, for example, a stress in-
tensity factor in fracture mechanics or a drag coefficient in fluid dynamics.
Other goals may be the determination of certain critical data, such as
collapse points or buckling points and their associated loads in structural
mechanics. Sometimes special techniques, e.g., variational techniques, can
yield highly accurate estimates of particular functionals.
In connection with most of these goals we are interested in quantita-
tive results that have a desired accuracy. For this the error has to be
defined, in that a family of exact results has to be specified with which
the computed data are to be compared, a norm has to be prescribed in
terms of which the error is to be measured, and some procedure has to be
established for estimating the error.
Such error estimation capabilities are certainly important hi many
applications, provided they can be guaranteed to be reliable. For example,
for the many types of certification computations required hi the design of
complex structures or nuclear plants, the availability of reliable estimates
of the accuracy of the computed data is obviously important. In other
cases, such estimates may reduce the total design effort and avoid un-
necessary overdesign. At the same time, the availability of effective er-
ror estimates introduces the possibility of applying adaptive techniques to
structure the computation to achieve a desired error tolerance at minimal
cost or to provide the best possible solution within an allowable cost range.
Many of the theoretical studies of solution algorithms for classes of
64
differential equations provide for some a priori error estimates, usually
of an asymptotic type. For the practical error assessment these a priori
bounds are rarely computable or very accurate. Thus, one is led to the
necessity of developing a posteriori error estimates that utilize information
gained during the computation itself.
For the finite-element solution of certain classes of elliptic boundary-
value problems, some computable and reliable a posteriori error estimates
have been developed and analyzed in recent years. Most of these a pos-
teriori estimates are based on a local analysis. For a given mesh A consist-
ing of elements AI, . . ., A n and with a corresponding approximate solution
on that mesh, an error indicator rjj is associated with the yth element A J(
These r\j have to be computable in terms of information about the prob-
lem and the approximate solution on Aj and, at most, on the immediate
neighbors of that element. On the basis of the indicators r}i,...,r) n an
error estimate e(A) is then constructed. Of course, the T)J as well as e(A)
depend on the chosen norm. The theoretically important question is then
the relationship between e(A) and the norm of the actual error e(A). The
effectivity of the estimation may be judged by the effectivity index 9 =
(A)/||e(A)||. In practice, it is usually more important for 6 to be close
to 1 than that < 1. Moreover, it is essential that 6 converges to 1 when
||e(A)|| tends to zero, so that for an accuracy of, say, 5 to 10 percent the
value of \0 1| is expected to be less than, say, 0.1 or 0.2.
The requirement of designing error estimators with these properties
for realistic classes of problems and various different norms certainly repre-
sents a demanding research task. The results available so far suggest that
such error estimators can be developed at least for linear problems. For
nonlinear problems the situation is in an embryonic stage although some
results available for model problems in one space dimension indicate that
estimators for the error along continuation paths and for the location of
critical points are computationally feasible. There is certainly considerable
need for concentrated research in the general area,
3.8.2 Other Validation Measures
Checking for mesh dependencies and convergence rates ie another way to
help validate a model. If a method is supposed to be, say, "second order,"
then as the mesh intervals are halved, the errors should be reduced by a
factor of 4. If this does not happen according to one's error estimate, then
the method does not have the desired order yet and may be in error. Also
one should continue to refine the mesh until the results are independent of
refinement to within the desired accuracy. Solutions can also be checked
65
for sensitivity to mesh orientation. The solution should not be aligning
itself with some special property of the mesh network.
Numerical methods can also be validated by comparing them with
analytic solutions. Often this is possible by choosing limiting cases where
an asymptotic analysis applies. In some cases complex constitutive laws
can be replaced with constants so an exact solution is possible. Higher-
dimensional models can sometimes be validated by comparison with pre-
viously confirmed lower-dimensional results. Bounds on solutions are
another way to help verify a model. Often conservation laws or minimum
and maximum principles apply, or the solution is known to approach a
previously known equilibrium or steady-state condition. These properties
can be used as a basis for validation.
To the extent possible one should always write computer codes in
a modular fashion such that each part of the model can be validated
separately. Even though this is just good structured programming prac-
tice, and not necessarily related to mathematics, it is an important aspect
of model validation.
3.8.3 Sensitivity Analysis
In addition to the independent and dependent variables, most models also
depend on certain physical parameters. Unfortunately, those parameters
are often not known accurately. Sensitivity analysis is a systematic
means to quantify the relationship of model parameters and model results.
Doing the sensitivity analysis requires solution of an additional system
of differential equations. These equations are a formal statement of
the relationship between the dependent variables of the system and the
parameters. The results are given in terms of a matrix of partial deriva-
tives of the dependent variables with respect to the system parameters.
One way to think of the analysis is as a method to provide theoretical
"error bars" for the model.
New methods are being developed to solve these sensitivity equations
quite efficiently. The methods rely on the observation that the sensitivity
equations are linear, regardless of the nonlinearities in the original model.
The method solves differential equations for the Green's function of the
sensitivity equations, and then the sensitivity coefficient matrices are com-
puted by quadrature for the various imhomogenous terms corresponding
to the parameters. These methods have been successfully applied to prob-
lems in chemical kinetics.
66
4. NUMERICAL METHODS
4.1 DISCRETIZATION METHODS
Many important physical problems are modeled by boundary or initial-
boundary-value problems. In this modeling, the physical state under
consideration is characterized by a function u, which is the unique solution
to the boundary or initial-boundary-value problem. Thus a major part of
the analysis of the physical problem is the determination of u or some of
its derivatives or some functional of it.
In nearly all important problems the determination of u is an infinite-
dimensional problem in the sense that u does not lie in an explicitly
known finite-dimensional space or, alternatively, cannot be expressed in
terms of a set of explicitly known functions by means of a finite set of
parameters. Thus u, cannot be "completely" represented on a computer,
and it is necessary to resort to some type of approximation or discretization
method. In essentially all discretization methods one attempts to construct
a function, u &pprox , which on the one hand is "close" to u and on the other
is characterized by a finite set of parameters, and so can be calculated.
The aim of finding ii approx that is "close" to u is made precise by
requiring that \\u - u &pprox \\ < r (or that \\u - u appr ox || < r||u||), where
1 1|| is a physically relevant norm and r is a physically relevant tolerance.
The goal is thus to find u app rox satisfying this error criterion with the
least expenditure of computational effort. The selection of a discretization
procedure is influenced by a number of considerations, the most important
of which we now list:
1. The goals of the analysis of the physical problem.
2. Known mathematical properties of the physical problem and the
algorithm.
3. Hardware considerations, e.g., the availability of parallel process-
ing.
4. Computer-science considerations, e.g., data-management require-
ments,
5. Restrictions on computation time and expense.
67
As indicated above, the exact solution is approximated by a function
u a pprox , which is expressible in terms of a set of explicitly known (basis)
functions by means of a finite set of parameters, which are determined
in the computation. One broad classification of methods is in terms
of the nature of the basis functions, namely (a) those involving global
basis functions having global support and (b) those involving local basis
functions having small support. Another classification is based on tne
extent to which the method is adaptive. Adaptive methods refine and
modify themselves on the basis of partially completed computations.
Of course, these classifications are not precise, and there are methods
possessing the different properties in various degrees. We now turn to a
more detailed discussion of the most commonly used discretization proce-
dures.
4.1.1 Finite Differences
One of the most frequently used discretization methods is the method
of finite differences. The central idea is to replace each partial deriva-
tive occurring in the differential equation in the underlying boundary or
initial-boundary-value problem by an approximating difference quotient,
For example, the first-order derivative du(x,y)/dx may be replaced by
the forward difference [u(x + h,y) -u(x,y)]/h or the backward difference
[u(x,y] - u(x - fr,y)]//i, whereas the Laplacian Au = u xx + Uyy may be
replaced by the five-point difference operator
u(x + h.y} + u(x. v + h) + u(x- h.
This replacement of all derivatives in the differential equation by ap-
propriate difference quotients leads to a system of equations, called
difference equations, for the numbers Uij, which are to approximate the
values of the exact solution u(x,y] at the finite-difference mesh points
(yhi kh), j, k = 0, 1, +2, . . ., where h is a small positive number, called the
mesh parameter. Nonuniform meshes may be considered as well, but at
considerable complication in the form of the difference quotients in the
resulting difference equations.
There are a number of important questions that arise in connection
with the study of difference methods, namely questions of
68
1. Accuracy of the method;
2. Stability of the method;
3. Efficient solvability of the finite-difference equations;
4. Robustness of the method with respect to the input data (e.g.,
coefficients, forcing functions, meshes); and
5. For nonlinear problems, the separation of legitimate approximate
solutions, which correspond to an exact solution, from spurious ap-
proximate solutions, which do not correspond to any exact solution.
Finite-difference methods are an example of the local basis class of
approximation methods. They can be adaptive. Adaptivity is usually
introduced via adaptive mesh selection: the mesh chosen at any stage in
the computation is based on the previous computations.
Much important work on difference methods remains to be done. We
mention in particular the construction of effective difference methods for
nonlinear problems, in particular for those problems whose solutions have
shocks.
4.1.2 Variational Methods of Discretization
We will now discuss a class of discretization methods based on variational
formulations of the physical problem under consideration. As noted above,
the problem of calculating the exact solution u is infinite dimensional in
the sense that u is only known a priori to lie in an infinite dimensional
space, say H. For a linear problem, a variational method of discretization
consists of (a) a finite-dimensional space S cH, called the trial space,
in which the approximate solution is sought, (b) a finite dimensional test
space V, and (c) a bilinear form B(u, v) defined on H X V. The approximate
solution is then determined by requiring that
u
approx
-B(ii a pprox > i>) = B(u, 11), for all v 6 V
where for veV , B(u, v) is computable from the data of the problem without
knowing u, Since and V are finite dimensional, u app rox can De calculated
by means of the solution of a system of linear equations if the original
equation is linear. Usually, for nonlinear problems, B is nonlinear in u,
and the resulting system of equations becomes nonlinear.
69
The approximate solution depends on the selection of S,V, and B.
For any problem there exists a wide variety of variational methods of
discretization, i.e., a wide variety of choices for S, V, and B. The rational
selection of S, V, and B is a central problem.
Variational methods can be of either the local-basis or global-basis
type. Of those of global-basis type we mention the various versions of the
spectral method and the p version of the finite-element method. Finite-
element methods, considered in more detail below, and collocation methods
are typically of the local-basis type. We also note that variational methods
can be adaptive.
4.1.3 Finite-Element Methods
Finite-element methods arise if for 3 and V we choose spaces S h and V h
of piecewise polynomials of fixed degree. That is to say, the underlying
spatiaMomain for the problem is broken up into small geometric pieces,
called "elements," whose size is measured by a parameter h. The functions
in both S h and V h are then restricted to be polynomials in x and y on
each piece but allowed to be different polynomials on the different pieces.
In designing finite-element methods, i.e., in selecting S ht V h , and B, one
attempts to achieve approzimability, stability, and systems of equations
that can be solved effectively.
Approximability here refers to the ability of the space S h to ap-
proximate the unknown solution u. The solution u is unknown a priori,
and often only the information u e H is available. In such situations
6 h has to be selected so that every function in H can be approximated
sufficiently well by one in S h . However, an approximation based on a few
large elements can provide additional information on u, which can be used
in turn to refine those elements. This type of adaptive element selection
is especiaUy important for problems with sharply varying solutions (see
section 4.2).
In the choice of the bilinear form B one has, in effect, the freedom to
choose a variety of variational principles, many of which have a natural
physical connection with the original problem. For instance, the stan-
dard Rita method is based on the principle of minimum potential energy.
An alternate variational principle is the complementary energy principle.
Approbation methods based on this principle are referred to as com-
plementary energy or equilibribum methods. These methods involve a
constraint that is usually difficult to satisfy. One way to circumvent this
difficulty is to treat the constraint by means of a Lagrange multiplier. This
70
leads to the so-called mixed methods. They appear to be promising for
many important problems and have recently received a large amount of
attention.
4.1.4 Transformation Methods
In the transformation or pseudospectral method, the discrete approxima-
tion u is first mapped by a transformation of the form
into the m-dimensional function space of the coefficients, a,j. The basis
functions fa and the transformation are chosen so that T and T~ l are
fast (order of m log m operations) and so that differentiation D is simple
in transform space. The derivative approximation can then be written as
dx
Some common transforms are based on the fast Fourier transform, where
the fa are trigonometric functions or Chebyshev or Legendre polynomials.
Selecting the <t> 3 - as piecewise polynomials with compact support, such as
the B-splines, is another good choice. By choosing the transformation
to incorporate some crucial property such as the periodicity or symmetry
of u one can improve the accuracy of the method for a fixed number of
basis functions. This can sometimes best be done by choosing the basis
functions close to the eigenfunctions of the differential equation.
4.1.5 Semidiscrete Methods
When solving a partial differential equation one sometimes discretizes with
respect to some but not all of the variables. For example, for the diffusion
equation governing the cooling of a hot rod, one may discretize with respect
to the space but not the tune variable, thereby replacing the original
partial differential equation by a system of ordinary differential equations.
Such an approximation method is referred to as a semidiscrete method or
as the method of lines. Semidiscrete methods may be used for hyperbolic
as well as parabolic equations.
71
One has, of course, eventually to discretize with respect to the time
variable as well, Semidiscrete methods are based essentially on the view
that very effective time discretization methods are available (in the form of
preprogrammed software packages for the solution of ordinary differential
equations) and that the spatial discretization is the main concern. An
alternative point of view is to consider both discretizations simultaneously.
Such fully discrete methods have been analyzed and tested extensively.
Among important research problems for semidiscrete methods we
mention the problem of adaptive mesh selection for the spatial discretiza-
tion.
4.1.6 Method of Characteristics
This is a method for hyperbolic equations, particularly for those involving
only one space dimension. In these equations the solution at some point
in space-time depends primarily on its values along a fixed, finite number
of curves (characteristics) going back in time from the given point. The
approximations to these values are determined from difference equations
that are closely related to the characteristics of the differential equations.
This method has a natural generalization to quasi-linear second-order
equations in two independent variables. It is especially important for
problems whose solutions have shocks.
4.2 ADAPTIVE GRID METHODS
For realistic problems it is rarely feasible to design numerical processes that
are rehable and accurate at a reasonable cost and yet that do not utilize
some form of adaptivity. Put simply, most two-dimensional and virtually
aJI three-dimensional problems are undercomputed without this technique.
This statement will almost certainly remain true after the next generation
of computers is available. The adaptive approach is to distribute the com-
putational effort nonuniformly, so as to concentrate on the most singular
parts of the solution or the most important parts of the solution from the
point of view of critical design parameters. Correspondingly one must
give less computational effort to the regular or less critical parts of the
SL T 'I** SeCti ^ 4 ' 9) - M ** 8ame time > >-P**VB approaches
may also lead to a simplification of the input data needed for the program
72
and hence free the user of part of the drudgery typical in preparing such
input and in having to make the many a priori decisions required by most
of today's programs.
The goal of adaptive grid methods is to choose a grid that is par-
ticularly refined, or that is aligned or oriented optimally, with respect to
the solution in regions of space and time that are critical for solution ac-
curacy. Thus adaptive grids utilize local mesh refinement or optimal local
mesh orientation.
The simplest adaptive grid is a preliminary, analytically determined
coordinate transformation. For example, in the transformed coordinates
the problem may be approximately independent of one variable or other-
wise simplified. The next strategy is to choose a numerically deter-
mined coordinate transformation. Typically, in two dimensions the coor-
dinate transformation is obtained by the solution of a subsidiary partial
differential equation. The resulting grids may be expected to give both
improved mesh refinement and mesh orientation. The method is some-
what problem dependent and can give rise to discretization errors in the
mapping of solution variables between the various grids.
Lagrangian grids for tune-dependent problems are adaptive for
material interface singularities, because the grid points move with the
material particles. Since these well-established methods also develop
rezoning and mesh entanglement problems, they provide a reservoir of ex-
perience concerning the difficulties associated with other evolving, adaptive
algorithms.
A refined grid can be constructed without recourse to a coordinate
transformation. In response to a critical solution parameter or solution
error criteria, selected regions of space can be flagged, preferred orientar
tions selected, and refined subgrids introduced locally. Then, small time
steps are chosen on the finer subgrid, and an interpolation problem must
be solved to blend fine and coarse grid solution values. Finally, the con-
struction is recursive, so that three, four, etc. levels of refinement can be
introduced automatically, in response to some a posteriori error estimate
on the next coarser grid. Precisely defined, reliable error estimators appear
to be central to the design of effective adaptive procedures.
In the context of finite elements for elliptic problems there is a devel-
oped theory for a posteriori error estimates (and, hence, for adaptivity),
which is based on local analysis. For a given mesh made up of elements
AI, . . ., A n and with the aid of the corresponding approximate solution, an
error indicator r)j is computed for each element Aj. For certain classes of
problems it has been proven that asymptotically the errors become optimal
when all the r\j become equal. This equilibration principle provides the
basis for the control of the adaptive process. In essence only those elements
73
A.,- are refined for which the error indicators are too large in comparison
with those of the other elements. And elements that are unnecessarily
small are combined with their neighbors. The study of suitable algorithms
for this has only recently begun. For example, some results have been ob-
tained for algorithms based on the assumption that during the refinement
process none of the elements could be combined without increasing the
maximum value of the error indicators.
In tune-dependent problems, differential equations can be introduced
for the ever-changing optimal location of the grid nodes. These equations
are then solved simultaneously with the original partial differential equa-
tion, leading to what are known as moving mesh methods.
Alignment of the grid without refinement is also possible, if fixed
or moving boundaries or interior interfaces are specified as part of the
problem. By alignment (without refinement), a regular grid index structure
can be preserved. Maintaining the grid structure has the advantage that it
potentially allows fast iterative methods, such as fast Fourier transforms,
to be used as part of the solution algorithm.
Adaptive grid techniques have been applied successfully to a range of
time-independent problems. The newer time-dependent methods need to
be developed to the point where they can be applied to meaningful two-
and especially three-dimensional problems and compared with alternative
methods. An important question requiring further attention is the con-
struction of efficient, computable a posteriori error estimates for realistic
classes of problems. Specifically, even for steady-state problems, there is
the question of the validity of the equilibration principle mentioned above
as well as the design and analysis of suitable adaptive control laws to
implement this principle. This latter problem may require examination
and incorporation of results in such fields as automatic control theory,
artificial intelligence, and learning processes. It also raises the problem
of the choice of appropriate refinement techniques that produce meshes
with desirable properties. Last but not least there is the question of data
management, which must be handled successfully to control the vastly in-
creased algorithmic complexity associated with adaptivity and to achieve
computational efficiency.
74
4.3 METHODS FOR SOLVING DISCRETIZED DIFFERENTIAL
EQUATIONS
Finite-difference and finite-element discretizations of partial differential
equations usually give rise to large systems of equations in which each
unknown is coupled to only a few of the other unknowns. Systems with
tens or hundreds of thousands of equations are relatively common. For
time-dependent problems these systems arise from the use of implicit
time discretizations. For sufficiently fine grids, the numerical solution of
these systems consumes a major part of the computer time for an entire
simulation.
Most nonlinear systems are solved by some form of iterative method
based on linearization, such as Newton's method. At each step, these
methods result in large sparse linear systems. In many cases, iterative
methods converge only if the initial guess is sufficiently close to the solu-
tion.
There are two basic approaches to solving large sparse linear systems
of equations: direct sparse matrix methods (i.e., some form of Gaussian
elimination that takes advantage of the location of most of the zeros)
and iterative methods, where the sparsity makes each iteration relatively
inexpensive. No method is best for all problems. For many problems
a combination of methods is attractive. Usually this takes the form of a
block iterative scheme using a direct method to solve the subsystems whose
diagonal blocks are matrices. Combination methods are of particular
importance for the solution of linear systems arising from coupled systems
of partial differential equations.
Direct methods are relatively well understood today, and a number
of excellent implementations of them exist for serial computers. Direct
methods vary in the extent to which they take zeros into account. The
simplest, nontrivial approach is band elimination, which takes account
only of those zeros occurring outside the band determined by the extreme
nonzeros of the matrix. The most complex approach is general sparse
elimination in which all the entries that remain zero during the alimrnflr
tion are taken into account. For systems that do not require pivoting
(e.g., those having symmetric, positive definite matrices) there is a sym-
bolic preprocessing phase in which the location of these zero entries is
calculated. There exist good techniques (the nested dissection ordering
and the minimum degree ordering) for arranging the unknowns and the
equations so as to minimize the zero fill-in during elimination. For model
problems it has been shown that the nested dissection ordering provides
asymptotically optimal results with respect to work and storage. The
75
minimum degree algorithm is a valuable heuristic approach that is com-
petitive in practice but that has defied analysis. For systems that require
pivoting for numerical stability one cannot compute the zero structure a
priori. Moreover, the ordering of the unknowns and equations to minimize
zero fill-in mil usually be significantly altered.
Some of the strong points of direct methods are as follows: (1) They
provide an exact answer (modulo round-off error) to the linear system with
a fixed number of operations independent of the system's condition num-
ber. Most production structural analysis packages such as NASTRAN use
some form of direct method, even for three-dimensional problems. In these
applications (many of which are based on fourth-order elliptic problems)
it is necessary to use higher-order precision because of conditioning prob-
lems. (2) For problems with two dependent spatial variables, their execu-
tion time is often shorter than the execution tune for iterative methods,
especially for problems of moderate size. Some of their principal disad-
vantages are the following: (1) they require considerably more storage than
iterative methods, even for two dimensions; (2) they will almost always be
noncompetitive with iterative methods for three-dimensional problems in
terms of running time; and (3) except for the simpler methods such as
band elimination (and to a lesser extent profile or envelope elimination)
they do not vectorize well. In fact, owing to the necessary indirect ad-
dressing involved, there are as yet no efficient implementations of general
sparse direct matrix methods for the CRAY-1 or CDC CYBER-205 super-
computers.
Except for structural problems, iterative methods are commonly used.
Classical iterative methods, such as successive overtaxation, Chebyshev
semi-iterative methods, and such newer methods as the preconditioned
conjugate gradient method, are fairly well understood for symmetric, posi-
tive definite systems, and they are easy to implement. However, the situa-
tion is not so bright for nonsymmetric or indefinite systems, though in
practice classical iterative methods may converge rapidly with a clever
choice of the iteration parameters. Nonsymmetric systems with indefinite
symmetric parts are especially difficult to solve, and no satisfactory general
algorithms are known at this time. Such problems arise in the simulation
of secondary and tertiary thermal recovery techniques for petroleum reser-
voirs (see Section 2.6).
Some of the principal advantages of iterative methods are the follow-
ing:
1. They tend to require minimal storage (proportional to the number
of unknowns);
76
2. They are reasonably fast for a wide range of problems. Moreover,
the number of iterations required is independent of the number of space
dimensions in the underlying partial differential equation (but not of their
order, or of the mesh size);
3. They can take advantage of good initial guesses to the solution, as
would be available in time-dependent or nonlinear problems; and
4. Many of them vectorize reasonably well on supercomputers.
Some of their disadvantages follow:
1. Mathematically rigorous stopping criteria may be difficult to for-
mulate, e.g., for linear systems with matrices that are not symmetric,
positive definite;
2. Many of the methods require a selection of iteration parameters,
and the performance of the methods depends critically on such a selection
(this difficulty is being overcome somewhat by the relatively new adaptive
methods and the preconditioned conjugate gradient-type methods);
3. The interaction between linearization and iteration is not well
understood, especially for discretizations of coupled systems of partial
differential equations; and
4. Nonsymmetric or indefinite problems may cause great difficulties
for iterative methods.
The relatively new multigrid iterative method combines the well-
understood behavior of a given iterative technique with the fact that an
underlying differential equation is being solved approximately. It alter-
nates iterating, toward the differential equation's solution, on fine and
coarser subgrids of the spatial network, with the goal of performing no
more computational work on the finer (hence expensive) grids than is ab-
solutely necessary. In many cases of practical interest, such as in neutron
diffusion in complex environments, the technique yields sufficiently ac-
curate solutions to the equations in a computational time proportional to
the number of unknowns. This has been a long-sought-for goal in the
approximation of elliptic equations in more than one space dimension. As
an iterative method it also has a natural extension to nonlinear equations;
and its logical structure, together with the already necessary calculation
of residual errors, points toward the incorporation of adaptivity for the
nesting spatial grids.
77
4.4 CONTINUATION AND HOMOTOPY METHODS
For more than a century the so-called continuation principle has proved
to be an important tool in mathematical theory. For example, early uses
date back at least 100 years to H. A. Schwarz and his work on con-formal
mappings; then, in the early part of this century, J. Hadamard and M
Levy applied these techniques in connection with the inversion of nonlinear
mappings, and it is also a basic tool in J. Leray and J. Schauder's work on
differential equations. But it was essentially only in the early 1950s, with
the advent of automatic computing, that continuation approaches began
to be used in numerical applications. Since then they have grown into an
extremely powerful technique for the numerical solution of wide classes of
nonlinear problems.
One of the problems to which continuation techniques are applied
concerns the solution of a nonlinear equation F(x) = in some space
X. In order to compute a specific solution x* X a possible approach is
to imbed the equation into a family of equations H (x, t) = 0, < t < I,
for which there exists a computable solution path x = x(t), 0<t<l, in
X beginning at a known point z(0) and ending at the desired solution
z(l) = x*. In other words, one considers a family of smoothly changing
problems, the final problem being the original problem in question and the
initial problem being one whose solution is easily determined. Use of the
continuation, therefore, requires an ability to solve a sequence of problems
when the solutions to nearby problems are known.
A related, but conceptually distinct, problem arises hi many applica-
tions in science and engineering where the equation under consideration
always includes a number of physically relevant, intrinsic parameters. In
other words, the equation has the generic form F(x,p) = 0, where x belongs
to some state space X and p varies in a parameter space P. In this setting
it is rarely of interest to determine only a few specific solutions. Instead,
the requirement is to follow paths on the solution manifold {(z,p) e X X
P} F(x, p) = 0} in the space XxP of all state and parameter variables and
thereby to detect specific features of the manifold that signify a change of
behavior of the system under study.
For the first of these problems, namely the computation of a specific
solution of some question, two distinct continuation techniques are avail-
able. The first involves the case when the path of solutions to the family
is smooth, which in turn allows its representation as a solution of a
differential equation. The second approach is based on homological rather
than homotopic concepts and makes use of a numerical form of a result in
algebraic topology (Speraer's lemma). This approach has been refonnu-
78
lated and is now considered principally in the form of algorithms involving
piecewise-linear solution paths. Much of the recent research in this area
concerns these piecewise-linear continuation methods and their application
to fixed-point problems in economics optimization and game theory.
Continuation methods for following paths on the solution-manifold of
parametrized equations developed mainly in structural engineering under
the name of incremental methods. Applications to other areas, as, for
example, to computational fluid dynamics and to phase transitions in
statistical physics have only begun to appear relatively recently. For a
numerical analysis of a given solution, manifold continuation methods have
to be considered hi a broader sense as a collection of numerical procedures
for completing a variety of tasks, including the following:
1. Follow numerically any curve on the manifold specified by an a
priori given combination of parameter values with one degree of freedom.
2. On any such curve determine the exact location of target points
where a given state variable has a specified value.
3. K desired, at any such target point switch to the trace of a different
solution path;
4. On any such curve identify and compute the initial points where
stability may be lost;
5. From any one of the critical points follow a path in the critical
boundary. In contrast to the case under 1, such paths are typically
specified by combinations of the parameters with two degrees of freedom
together with the implicit requirement that all points of the path are in
the critical boundary.
6. On any solution path determine the location of bifurcation points
and the paths intersecting at that point.
For specific applications additional tasks may arise. For instance, if
the parametrized equation represents the equilibrium equation of a system
of differential equations then we may wish to locate Hopf bifurcation points
on a particular solution path where periodic solutions of the dynamical
system branch off from the static equilibrium.
In recent years much research has been devoted to these various tasks,
but there remain many open questions especially hi connection with the
more specialized tasks, such as the location of Hopf bifurcation points.
Moreover, library-quality packages for the basic tasks 1-3 are still under
development, and the software is by far not so advanced as, say, in the case
of software for the solution of ordinary differential equations. For specific
79
applications, as, for example, to fluid-dynamic problems or to the case
of structural problems involving plasticity, creep, and viscoelastic effects,
the situation is still wide open; and only relatively ad hoc techniques are
available.
4.5 OPTIMIZATION METHODS
Optimization problems occur in all areas of science, engineering,
economics, and applied statistics. They may involve some least-squares
approximation of observed data, fitting of parameters occurring in a math-
ematical model on the basis of experimental observations, optimization of
the design of an engineering structure, optimal control of an engineering
system, or modeling of economic or business systems. These are only a
few examples that do not even begin to cover the numerous types of op-
timization problems that arise in practice.
Broadly speaking, in all of these problems a real function, usually
called the objective function, is to be minimized or maximized over some
constraint set in a given finite- or infinite-dimensional space. The problems
differ considerably depending on the mathematical characteristics of the
objective function and the constraints, the dimension of the space, the
amount of computable information that is available, and the requirements
of specific applications.
A good overview of the research problems concerning computa-
tional methods for such optimization problems in finite-dimensional spaces
has been given in the report Program Directions for Computational
Mathematics, June 1979, prepared for the Applied Mathematical Sciences
Research Program of the U.S. Department of Energy (edited by R. E.
Huddleston, Sandia National Laboratories). This report identified the fol-
lowing research areas:
1. Analysis and comparison of techniques for dealing with nonlinear
inequality constraints.
2. Construction of algorithms for environments where computer
storage is restricted.
3. Production of high-quality software and other software-directed
activities.
4. Development of algorithms and software for problems in which
some of the variables are discrete or integer.
80
5. Techniques to assist in finding global (or specified) minima.
6. Large-scale linear programs with special structure, e.g., staircase
and block-angular.
7. Methods for nondifferentiable optimization.
8. Techniques for selected nonlinear estimation problems, such as
separable or constrained least squares, and data fitting in special norms.
9. Investigation of the connection between fixed-point methods and
alternative methods for solving nonlinear equations.
10. Further study of selected aspects of linear constrained problems,
such as special types of constraints, application of conjugate gradient
techniques, the merits of various active set strategies, and the influence
or scaling.
The report does not address infinite-dimensional problems, for ex-
ample, problems in control theory or the calculus of variations. Here the
mathematical problems often become formidable, and special techniques
are needed in most applications. For the computational solution, the prin-
cipal approach applies finite-dimensional optimization methods to a dis-
cretized form of the problem. This, in turn, raises the usual questions
concerning the convergence and quality of the approximations.
4.6 GRAPH-THEORETICAL METHODS
Many physical systems, or their mathematical models, involve a number
of discrete objects that are the producers or users of certain data or
commodities and that in turn are interconnected by links that can transmit
these items from one of the objects to another one. The underlying
mathematical structure is then a graph consisting of a collection of vertices
connected by arcs or branches.
The modeling of some classes of physical systems in terms of graphs
is rather natural, for example, in the case of communication or transpor-
tation networks where the branches represent, say, power lines, roads, air-
line routes, or water pipes, and, correspondingly, the vertices are power-
generation stations, communities, airline terminals, or water reservoirs. In
other cases the connection between a particular problem and its graph-
theoretical interpretation is less obvious, for example, in the case of a
discretization of a particular boundary-value problem on a given grid or
the representation of a collection of data in a computer.
81
Corresponding to the wide range of applications, the forms of graph
problems differ widely. For certain network models interest may center
on the connectivity properties of the graph, that is, on the determination
of whether a particular commodity can be transported between specific
vertices. Then questions arise about the maximum possible flow that can
be accommodated under particular constraints. These problems in turn are
related to the so-called vulnerability and reliability problems for networks.
On the other hand, if time enters into consideration then questions of
waiting time and of best routing may arise. These are only a few of many
types of such problems.
From a computational viewpoint these various problems have stimu-
lated the development of numerous combinatorial and graph-theoretical al-
gorithms. But there remain many open research questions, especially when
it comes to the production of general, robust software and the availability
of algorithms for problems involving large graphs.
Graph-theoretical formulations are also finding increasing application
in connection with the numerical solution of problems that do not have an
inherent graph-structure. One such class of problems concerns the com-
putational handling of large sparse matrices. Many of the algorithms used
here perform a sequence of steps that transform the matrix into succes-
sively simpler matrices of, say, a more nearly upper triangular or diagonal
type. These transformations achieve their aim by introducing zeros in
place of originally nonzero matrix elements. But, as an unavoidable by-
product they also introduce nonzero elements hi places where the original
matrix elements were zero. Thus, in order to exploit the sparsity struc-
ture of the matrix one must provide a data structure that either allocates
from the outset space for all the fill-in during the computation or that
allows for a dynamical allocation of space for the fill-in when it occurs. In
both cases, graph-theoretical approaches form the basis for the design of
desirable algorithms.
An example of the use of a static data structure is the case when a row
and column representation can bring the matrix into banded form with
small bandwidth. This is frequently the case of the linear systems arising
hi the finite-element discretization of elliptic boundary-value problems,
and the corresponding bandwidth optimization routines are widely used
in practice. On the other hand, if a dynamic storage structure is used, then
special care has to be taken in the design of pivoting strategies for reducing
fill-in while at the same time maintainmg numerical stability. Various
sparse matrix packages have been developed for this purpose. They are
typcially rather large and complex pieces of software. But the field is still
under active development, and there remain many open research questions.
82
Iii particular, much still needs to be done on the interrelation between
sparse matrix algorithms and special-purpose computer architectures.
In recent years graph-theoretical approaches have also found increas-
ing application in the design of algorithms for the numerical solution of
problems of mathematical physics, in particular, of fluid dynamics. For ex-
ample, it has been observed that certain natural finite-difference discretiza-
tions of the equations of viscous, incompressible flow admit interpretations
as systems defining flows on certain associated networks. Typically, the
network nodes correspond to the idealized control volumes represented by
the mesh points of the finite-difference equations, whereas the network
arcs correspond to the paths on which one may identify the discrete finite-
difference mass fluxes. Such observations can lead to remarkable savings
in computing costs because they open the way to a priori transformations
of the original (discretized) system of equations to completely equivalent
systems in substantially fewer variables. For flow problems in two and
three space dimensions the reduction factors are nominally 1/3 and 1/2,
respectively. In graph-theoretical terms this approach corresponds to the
transformation of the Maxwell-node equations to the Maxwell-mesh equa-
tions, long used in electrical-circuit problems.
The key requirement hi the application of network techniques to fluid-
flow problems is that the discretized momentum and continuity equations
be interpretable as "Ohm's laws" and "Kirchhoff-node laws" on an as-
sociated network. This requirement permits great generality in the actual
form of these laws and easily accommodates both implicit time-dependent
as well as nonlinear steady-state difference forms of the Navier-Stokes
equations.
It is even possible to extend the network approach to compressible flow
problems. Here the new idea appears to be the introduction of pseudo-
flows into the node law and the identification of corresponding pseudo-
characteristics to augment the Ohm's laws obtained from the momentum
equations. In the nonstationary case, this produces numerical schemes
whose stability is independent of local acoustic velocities.
The entire area is still under active research development. For ex-
ample, there are efforts to bring multispecies flow problems under the
purview of network theory. This includes, as an important special case,
two-phase flow problems that are a central concern in the design of reliable
energy-generating systems and the safety analysis of nuclear power plants.
83
4.7 SPLITTING METHODS AND DEFECT CORRECTIONS
Splitting is a means to separate a large intractable problem into a sequence
of smaller, or at least more easily solved, problems. These methods are
invoked to reduce significantly the computational effort (time and memory)
needed to solve a problem compared with solving it directly. Often a
form of splitting is required to make solution practical, and some cases
even tractable. Various forms of splitting are common in engineering and
scientific applications, even if they are not always recognized as such.
Defect correction is a widely used, if often unlabeled technique of
splitting. It presumes that one wants to solve a given hard problem, that
one has in hand a guess at its solution, and that one also has a nearby
problem that can be solved easily. It corrects the guess by solving the
easier problem with special data computed from the guess.
To illustrate the method, suppose that after discretization it is neces-
sary to solve a finite-dimensional system of equations of the form
-
where A is a nonlinear operator, 6 is a known vector, and v is the solution.
Often v is difficult to obtain directly, but the residual error
r A(w) b
for an approximate solution w is easy to evaluate. If there is a related
system
that approximates the original system and that is easier to solve, the defect
correction algorithm may be appropriate. Given an estimate v n near a
root v n+ i of the original system, we can expand the original equation in
a Taylor series to get
AK) - 6 + P(v n+1 ) - P(v n ] - ( J P - J A )(v n+l - v n ] + 0(e 2 )
where e = v n+ i v n . The defect correction iteration is any approximation
to the above equation. The simplest such iteration is
This iteration will converge if v n and J P (the Jacobian of P) are near
enough to v n+i and J At respectively.
84
One of the most common splitting algorithms in engineering applica-
tions is the alternating direction implicit (ADI) method. Here, the model
problem is a partial differential equation system having two or more in-
dependent spatial coordinates. Its direct solution requires too much com-
puter storage and time to be tractable. Using an ADI splitting, the prob-
lem is reduced to a manageable sequence of one-dimensional problems.
In terms of the defect corrections, the system P is the sequence of one-
dimensional problems that is much easier to solve, yet that closely ap-
proximates the original system.
Other forms of splitting also arise. For example, in certain systems of
partial differential equations some of the equations are weakly coupled to
the others. In these cases, solving the equations sequentially (rather than
coupled) can result in significant savings. Similarly, in some combustion
problems, considerable savings are realized through operator splitting al-
gorithms in which the chemical rate terms are handled separately from the
species transport terms. These methods are equivalent to matrix splittings
of the Jacobian of the system.
Analysis of splitting methods is important because the methods often
do not converge. We must be concerned about accuracy and convergence
of the factored system. Even though each iteration may be fast, we
must have some idea about the overall convergence rates and about any
degradation in accuracy resulting from the split. Analysis will also likely
lead to the identification of matrix properties that suggest a beneficial
splitting that would not be apparent from physical reasoning. Splitting
can be considered as either a splitting of the original equations or as an
approximate factorization of the iteration matrix. The former is the most
intuitive, and the one in which physical insight is valuable. However, the
latter is probably the one more amenable to analysis.
4.8 MONTE CARLO TECHNIQUES
Mathematical solution methods can be broken into two broad classes,
Monte Carlo methods and deterministic methods, depending on whether
chance events are included in the method. All scientists are familiar
with deterministic methods but most have little familiarity with Monte
Carlo methods. Deterministic methods do not depend on chance, and
calculations performed using the same input data and the same computer
code will provide exactly the same results. Monte Carlo calculations using
the same input data and the same computer code will provide different
85
results, depending on what chance events occur. In this section we discuss
what Monte Carlo methods are and how they might be improved.
The Monte Carlo method estimates averages on a probability modd.
Any quantity that can be expressed as an average value can be calculated
by Monte Carlo techniques. Sometimes a probability model is obvious
from the problem itself. For example, the probability that tossing three
fair coins and obtaining a desired outcome [for example two heads (H) and
one tail (T)] is easily estimated by Monte Carlo methods. The probability
model consists of assigning a probability of ^ to an H and a probability
of to a T on each toss and assigning a score of 1 to a desired outcome
(HHT, HTH, THH) and a score of to any other outcome (HHH, TTT,
THT, HTT, TTH). The computer has a random number generator that
generates random numbers uniformly on the unit interval (0,1). A uniform
distribution means that any random number is equally likely to have any
value between and 1. Thus a com toss can be simulated on a computer
by:
H occurs if < 5
T occurs if > f
To toss three coins the computer uses three random numbers &, &,
and 3 . A typical set of tosses might be i = 0.7 (T), & = 0.1 (H), and
3 = 0.4 (H), scoring 1. The probability of obtaining two heads and one
tail is (approximately) the average score (|) after many sets of computer
tosses.
Sometimes a probability model is not immediately apparent, but after
a little thought the desired calculation can often be expressed as the
estimation of an average value. For example, the integral
can be thought of as the average value of f(x) on the interval (a, b). To
see this, note that by definition the average value of f(x) with respect to
a probability density p(x] is
f b
</>=/ f(x)p(x)dx.
Ja
Thus J is the average value of f(x) with respect to the probability density
p(x} = l/(b a). Here p(x) is a uniform probability density on (a, b), and
86
x values can be samples uniformly on (a, 6) by setting
whereupon
i N
where I is the estimate of I. . , , . ,
There are many applications for the Monte Carlo technique, but the
application to neutron and photon transport has probably consumed more
computer time than all other applications combined. Because ol tms, ana
because neutron and photon transport problems have a natural pro baouiy
model, the remaining discussion will pertain to neutron or photon ^trans-
port problems. Furthermore, because neutron and photon transport are
similar, the remaining discussion will refer to particle transport.
Particle flow in nature is a random process. One cannot say exactly
what will happen to an individual particle. One can only say what the
probabilities are that various events occur. For example, a particle ttavei-
ing through matter has a certain probability to collide per unit s ,
the actual distance between collisions is unknown, but the probab iihty ol
traveling a distance I without collisions is known. Smulariy the ^nucuae
a particle will collide with is not known, but the probability of colliding
with the nuclide is known. ^,Wo-mQ ns
The simplest Monte Carlo model for particle transport P^lems uj
the natural probabilities that various events occur for the probabib
model in essentially the same way as the coin toss example. That
particles are followed from event to event in a ComputerLand tne
is always selected (using the random number generator) trom a
possible next events according to the natural event ******
type of model is called the analog Monte Carlo model because it is
analogous to the naturaUy occurring ^transport
particles contribute, and the statistical uncertainty m the answer
models for particle transport that will estimate the same average value
as the analog Monte Carlo model. These other techniques are known
as nonanalog Monte Carlo models, and they are useful because although
the average value remains unchanged, the variance (uncertainty) of the
estimate can often be made much smaller than the variance for the analog
estimate. Practically this means that problems that would be impossible
to solve in days of computer time can be solved in minutes of computer
time.
A nonanalog Monte Carlo model attempts to follow "interesting" par-
ticles more often than noninteresting particles. An "interesting" particle,
by definition, is a particle that contributes a large amount to the quantity
(or quantities) that needs to be estimated. There are many nonanalog
techniques, and they all are meant to increase the odds that a particle
scores (contributes). To ensure that the average score is the same in the
nonanalog game as in the analog game, the score is modified to remove
the effect of biasing (changing) the natural odds. Thus, if a particle is
artificially made q times as likely to execute a given random walk, then
the particle's score is weighted by (that is, multiplied by) 1/q. The average
score is thus preserved because the average score is the sum, over all ran-
dom walks, of the probability of a random walk and the score due to that
random walk. The trick in obtaining low- variance solutions is to choose g's
such that all random walks contribute the same score, in fact the average
score. This then would be a zero-variance solution.
It is always possible for any (linear) particle transport problem to
select q's for each random walk such that every particle's score is the
average score; that is, a zero-variance solution. The hooker is, of course,
that the perfect q'& are not known, thus a zero-variance solution is not
practical. However, people have often been successful enough in guessing
good g's, that is biasing the odds, to improve the calculational efficiency
by several orders of magnitude. This is typically done iteratively with a
person making several short Monte Carlo computer runs (calculations) and
using information gained on each run to better guess how to bias the next
run. When the person's guesses no longer improve the calculation, a long
run is done with the optimal biasing learned in the short runs.
Can the computer learn to optimize the biasing? The computer is,
after all, supplying the information that the person uses to learn. In
the past decade a number of computer learning techniques have been
tested. Thus far it has proven impossible for a computer to take an
arbitrary transport problem and, without human aid, optimize the biasing.
However, two things are worthy of note. First, computer learning aided
by human judgment appears to be substantially better in many cases than
human learning alone. This typically works by having the computer inform
88
the person how the computer would bias the subsequent run and having
the person selectively alter the computer's suggestions according to the
person's intuition. Second, the amount of human judgment required is
decreasing. The day may come when the computer needs no human aid.
Once human aid is no longer needed, the computer can learn to
adjust its own biasing particle, as the calculation proceeds. An interesting
implication of an adaptive Monte Carlo technique is that the common
central limit theorem of statistics that would constrain the accuracy of
the calculation to decrease as the square root of the number of particles
followed no longer applies. The common central limit, theorem applies
only when each particle's random walk is independent of all others and
the sampling process is identical for all particles. Consequently, the rate
of convergence may be more rapid. Although some limited statistical
theories exist for dependent random variables, it has not been investigated
whether these theories can profitably be applied to Monte Carlo particle
transport problems. Thus, it is uncertain how fast an adaptive Monte.
Carlo calculation is converging or indeed what the maximum convergence
rate for a good learning process might be. Empirical evidence shows that
the convergence can be substantially faster than the square root of the
number of particles. In light of this, Monte Carlo methods can be expected
to become more competitive with deterministic calculations.
4.9 PROBLEM-DEPENDENT METHODS
A variety of adaptive methods have a common goal: to overcome the
computer size and grid resolution limitations, which are especially severe
in singular problems, by use of computational elements that fit or model
the singularity more directly. In this way, it may be possible to incorporate
into the solution algorithm considerable analytic information beyond that
provided by the equations themselves. This theme occurs in many aspects
of numerical analysis. In the finite-element method, one may choose
elements that include any singularities in the solution being computed (see
Section 3.3). Grid adaptivity (see Section 4.2) is also problem dependent.
Here we discuss examples of somewhat more special methods. Of necessity
they apply only to classes of problems that possess related singularities.
89
4.9.1 The Riemann Problem and Nonlinear Wave Modes
The Riemann problem for a nonlinear hyperbolic system of conservation
laws is the Cauchy problem in one space dimension for data that are
everywhere constant except for a single jump discontinuity at the origin.
The solution of the Riemann problem provides a resolution of this discon-
tinuity into the nonlinear modes, or waves, which propagate coherently
in time. This idealized problem can be thought of as an approximate
description of a higher-dimensional flow field especially in the neighbor-
hood of a discontinuity hypersurface. This point of view has led to a num-
ber of improved finite-difference schemes, which attempt to determine the
various nonlinear wave modes at each point in space and time and to ad-
just the differencing of the differential equations to take advantage of this
knowledge. This adaptivity is to the structure of state space, in contrast
to the coordinate space adaptivity discussed in Section 4.2.
4.9.2 Front Tracking
Front tracking is a combination of adaptive grid methods with the use of
Riemann problems. The method is adapted to problems that have impor-
tant singularity hypersurfaces (lines in two space dimensions, surfaces in
three space dimensions), such as shock waves and contact discontinuities.
In one version of this method, there is an overall time-dependent coordinate
transformation to map the discontinuity into a fixed grid line in some set
of computational coordinates. In another version, the discontinuity is rep-
resented by a moving lower-dimensional grid that follows ("tracks") the
dynamical motion of the discontinuity. The motion of the discontinuity
and of the moving grid points that track it are governed by solutions of
Riemann problems, or equivalently by the method of characteristics.
4.9.3 Vortex Method
This method introduces small lines or surfaces of vorticity into a fluid flow.
The method is adapted to the study of shear-layer discontinuities, bound-
ary layers, boundary-layer separation, and turbulence. The equations of
motion of an ideal fluid yield simple equations for the motion of a collec-
tion of vortices imbedded in the fluid. In fact, the vortices move passively
90
with the fluid, and their mutual interaction is described by a Hamiltonian
system of ordinary differential equations, with Coulomb type interaction
energy. La the case of the Navier-Stokes equations, the vortex motion also
contains a diffusion term.
These methods have been applied successfully to the problem of tur-
bulent flame propagation (see Section 2.2). In this problem, the turbulent
mixing is a primary factor in determining the flame speed. The turbulence
comes from the boundary layer and in the boundary layer is calculated by
a vortex method. Related methods have been developed under a variety of
names such as boundary integral methods and Green's function methods.
4.9.4 Scale Invariance and the Renormalization Group
Scale transformations are the transformations x * sx, acting on space or
on space and tune. A function u is homogeneous of degree a if
u(sx, ay, sz) = s a u(x, y, z)
and scale invariant if a = 0. Many problems have solutions that are scale
invariant or approximately scale invariant over some parameter range.
Such solutions are called similarity solutions. Using the scale invariance,
one of the independent variables can be eliminated, making the solution
more elementary to compute.
However, scale invariance can also indicate the occurrence of com-
plex phenomena. Specifically, any singularity that occurs in a scale-
invariant problem must be repeated in all length scales (for which the
scale invariance holds). Mathematically, Cantor sets, snowflake curves,
and fractals are examples of such phenomena. In statistical physics, criti-
cal phenomena in the equation of state is a scale-invariant problem. One
general picture of turbulence holds that scale invariance (vortices on a
large range of length scales) describe the inertia! range, or energy transport
range, of turbulence.
To implement these ideas in a computational algorithm, one integrates
over a given set of length scales in a standard manner. The result of this
computation is taken as data for a new calculation over a new set of length
scales (with the original degrees of freedom eliminated from the prob-
lem). This process is iterated and if convergent gives the scale exponent
a, Sample numerical calculations of this type were discussed in Section
2.1 in connection with turbulence. The method is well established for per-
turbative calculations (with a small parameter) in critical phenomena in
statistical physics.
91
25
10
30
15
35
20
37.5
FIGURE 4.1 The evolution and merger of isolated vortex structures as
predicted by contour dynamical techniques. [From E. A. Overman n and
N. J. Zabusky, Phys, Fluids 25, 1297-1305 (1982).]
92
Since the renormalization group methods are novel, it is worth men-
tioning that the mathematical foundations of this method have been estab-
lished recently in several cases including examples of hierarchical models
in statistical mechanics, interval maps, and renormalized quantum fields.
4.9.5 Contour Dynamical Methods
Contour dynamical methods are being applied to a variety of incompress-
ible flows in two dimensions. These generalizations of the "waterbag"
method provide simplified models for following the evolution of contours
separating regions of piecewise-constant (pc) density that are the sources
of the flow. The flow is the result of the self and mutual interaction of con-
tours that evolve, mainly by area-preserving maps. These methods have
been applied to the Euler equations, where pc finite-area-vortex regions
and/or vortex sheets at density interfaces are sources of the flow; and
the equations for a weakly ionized and strongly magnetized ionospheric
plasma cloud, where pc ion-density regions are sources of the flow. For the
former, a large class of steady-state translating and rotating solutions with
pc vorticity ("V-states") have been found. Figure 4.1 shows the merger
and breaking for a perturbed corotating V-state (two pc finite-area vortex
regions with the same circulation), a familiar process in free-shear layer
experiments. Notice how the two regions merge to form one region (by
snipping out the common boundary at t = 10) and then stabilize by eject-
ing vorticity in thin filaments. With these methods it is possible to obtain
detailed information about the regions because the dimension is reduced
by one. The curvature of the contour provides a predictive signature of
the evolving small-scale structure, e.g., the roll up of vorticity filaments,
etc.
4.10 COMPUTER SOFTWARE AND HARDWARE
As previously discussed, large-scale scientific calculations that tax the
resources of the most powerful computers will continue to be essential to
modern research and development efforts. To obtain long-term reliability
and stability of future applications codes, implementing and testing of
high-level numerical software should be coordinated. This will require a
93
strong research and development effort with cooperation supported among
applications programmers, the theoretical numerical methods researchers,
and the computer-science software developers.
Repetition of expensive, error-prone, and time-consuming coding of
commonly used methods should be avoided. Much of the current scientific
software now being developed is redundant. If the common elements of the
existing codes were available as modules, future applications programmers
could use these routines and eliminate much of their efforts. New software
is most readily accepted if it is compatible with existing techniques and
simple enough so that potential users can observe tangibly better results
in a trial run than those existing methods can produce. If such high-
level routines were available, they could perform many of the common
procedures found in applications codes, including grid generation, rezon-
ing, numerical interpolation, differentiation, and integration; they could
approximate differential boundary conditions and solve large, sparse non-
linear systems of equations.
An important goal is the machine independence of applications pro-
grams. This is a difficult task because methods tailored specifically for a
particular machine architecture will probably become more the norm than
the exception. We can, by keeping machine-dependent codes in libraries of
high-level software with standard user interfaces, strive to keep the user's
scientific applications codes portable. The underlying library routines will
be, necessarily, less portable because the architecture of the new machines
will certainly be different from that of today's supercomputers. To uti-
lize the inherent powers of the new machines we will have to re-examine
traditional methods and identify the better ones for a particular machine
architecture.
The continuing revolution in microelectronics is having a profound
impact on scientific computing. Indeed, it is likely to change our concept
dramatically of what scientific computations can and cannot be effectively
performed. Most certainly the impact of this revolution will be much
greater than, say, the impact that floating-point hardware had. Moreover,
while the costs of individual tasks will be greatly reduced, the domain of
scientific computations will be greatly expanded and "frontier computa-
tions" will continue to be expensive.
These changes are being brought about by a number of factors:
Individual components are becoming increasingly faster and smaller. Very-
large-scale integration of such components is not only reducing the size of
the packaged systems but also providing opportunities for customized, in-
formation processors. Also, the decrease in the cost and size of computer
memory implies that we can look for much larger memory systems. This
will obviate many of the existing difficulties with secondary storage.
94
Despite the fact that components are becoming faster, the limits of
raw machine speed are not visible, and further gains will eventually have
to be made by using clever architectures and algorithms. Some kind of
parallelism seems to be unavoidable. The programming issues involved in
using parallel machines are still not satisfactorily resolved. The automatic
detection of parallelism and the resulting scheduling of multiple resources
are important open problems.
Architectures, such as systolic arrays, based on specific subtasks can
increase the performance of systems involving these subtasks by several
orders of magnitude and clearly have a bright future. Algorithms contain-
ing compute-intensive subtasks that can be vectorized in this fashion have
a bright future. Because these architectures are hi general regular, the
algorithms that can be vectorized for such machines tend to be regular,
i.e., simple, nonadaptive, uniform-grid, low-order algorithms. It is clear
that there are nicely behaved problems for which these regular algorithms
on specialized machines will require significantly less time than algorithms
requiring fewer operations on serial machines. It is also clear that no
matter how fast the specialized machines are, there are problems that
are sufficiently difficult that more sophisticated algorithms are needed for
more general-purpose computers.
In order to bring about these advances in architecture, it is necessary
to involve practitioners of scientific computation in the design process.
Luckily, modem design automation tools should make it possible for in-
terdisciplinary design teams to successfully synthesize innovative special-
purpose systems, and automated fabrication facilities should make it pos-
sible for such systems to be built, debugged, and used.
These advances in machine architecture should also have a dramatic
effect on the design and analysis of algorithms for scientific computing.
Traditionally, such analysis has been based on (asymptotic) estimates of
the number of arithmetic operations. However, with these new architec-
tures it is quite likely that the running time of an algorithm will be more
dependent on the movement of data than on the number of arithmetic
operations. Thus, we need to develop new analytic models of complexity of
scientific algorithms so that such models give us useful information about
the relative performance of algorithms.
95
Appendix A
Persons Contributing to this Report
Ivo Babuska, University of Maryland
K. Binder, KFA Julich GmbH, Institut Festkorperforschung
Norman Bleistein, University of Denver
Thomas E. Booth, Los Alamos National Laboratory
Billy L. Buzbee, Los Alamos National Laboratory
William S. Childress, New York University
Alexandre J. Chorin, University of California, Berkeley
Lawrence D. Cloutman, Los Alamos National Laboratory
James P. Corones, Ames Laboratory, Iowa State University
Joel E. Dendy, Los Alamos National Laboratory
Paul P. Eggermont, University of Delaware
Alan J. Faller, University of Maryland
Paul R. Garabedian, New York University
John M. Greene, General Atomics Company
Eugene Isaacson, New York University
Jasper A. Jackson, Los Alamos National Laboratory
Malvin H. Kalos, New York University
Herbert B. Keller, California Institute of Technology
Heins-Otto Kreiss, California Institute of Technology
Edward W. Larsen, Los Alamos National Laboratory
Peter D. Lax, New York University
Joel Lebowitz, Rutgers University
Dennis R. Liles, Los Alamos National Laboratory
C. C. Lin, Massachusetts Institute of Technology
Oscar Manley, U. S. Department of Energy
Steven B. Margolis, Sandia National Laboratories
Jerrold E. Marsden, University of California, Berkeley
Oliver A. McBryan, New York University
Keith Miller, University of California, Berkeley
Cathleen S, Morawetz, New York University
M. Zahair Nashed, University of Delaware
J. Tinsley Oden, University of Texas
Edward A. Overman H, University of Pittsburgh
Edward Ott, University of Maryland
96
George C. Papanicolaou, New York University
Charles S. Peskin, New York University
Linda R. Petzold, Sandia National Laboratories
Thomas F. Porsching, University of Pittsburgh
Garry H. Rodrigue, Lawrence Livermore National Laboratory
Jacob T. Schwartz, New York University
Philip E. Seiden, IBM Corporation
David H. Sharp, Lawrence Livermore National Laboratory
Lawrence A. Shepp, Bell Laboratories
Mitchell D. Smooke, Sandia National Laboratories
Blair K. Swartz, Los Alamos National Laboratory
Yvain M. Treve, La Jolla Institute
Michael Steenstrup Vogelius, University of Maryland
Burton Wendrqff, Los Alamos National Laboratory
Norman J. Zabusky, University of Pittsburgh
97
511.3 N27c c.l
National Research Council
(U.S.).
Computational modeling and
mathematics applied to the
University Libraries
Carnegie- Mellon University
Pittsburgh, Pennsylvania 15213
DEMCO
Live Music
Archive
Librivox
Free Audio
Metropolitan Museum
Cleveland
Museum of Art
Internet
Arcade
Console Living Room
Open Library
American
Libraries
TV News
Understanding
9/11